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Discrete-time Signals and Signal Processing 


A World of Signals and Signal Processing 


Technological innovations have revolutionized the 
way we view and interact with the world around us. 
Editing a photo, re-mixing a song, automatically 
measuring and adjusting chemical concentrations in 
a tank...each of these tasks requires real-world data 
to be captured by a computer and then manipulated 
digitally to extract the salient information. Have you 
ever wondered how signals from the physical world 
are sampled, stored, and processed without losing 
the information required to make predictions and 
extract meaning from the data? Signal processing 
is the study of signals and systems that extract 
information from the world around us. 


Signals, Defined 


Perhaps the best place to start the study of signal 
processing is with a dictionary definition: 


Signal* 


A signal is "is a detectable physical 
quantity...by which messages or information 
can be transmitted." 


*"Signal." Merriam-Webster.com 


This information aspect is critical to us. In other 
words, signals carry information. Signals are all 
around us; we encounter them throughout our day. 
Speech signals, for example, transmit language from 
one person to another via acoustic waves. If you're 
interested in looking for airplanes or other targets 
and sensing them by electromagnetic waves, you 
can use radar signal processing. Electrophysiological 
signals carry information about processes that are 
going on inside our bodies; tools like EKGs or MRIs 
work with those signals. And finally, financial 
signals transmit information about events in the 
economy, signals like stock prices over time or other 
economic markers. 


Signals are Functions 


A signal is a function that maps an independent 
variable to a dependent variable. Mathematically, 
we're going to think of signals as functions, and a 
function is simply a mapping from an independent 
variable that we can change onto a dependent 
variable that depends on that independent variable. 
The terminology x[n] is how we're going to denote a 
signal. It consists of an independent variable, n, that 
for each of its values produces the another 

value x[n]. 


Discrete-Time Signals 


Perhaps you are wondering about the use of 
brackets--instead of parentheses--in our signal 
function notation x[n]. This is the typical way to 
refer to discrete-time signals. A discrete-time signal 
is a signal where the independent variable n is an 
integer (as opposed to a continuous-time signal x(t), 
whose independent variable t is a real number). 


Plotting Discrete-Time Signals 


Recall that a discrete-time signal is a function 
with an integer-valued independent variable n. The 
variable n marches through time from negative 
infinity to positive infinity. For each value of n, we 
get the value of our function x[n]. Now, that x[n] 
is either going to be a real number, meaning it's 
going to live in the real number set, or it's going to 
be a complex number and live in the complex 
number set. 

Example of a discrete-time signal. 


Plotting Real Signals 


We're going to see a lot graphs in our study of signal 


processing: 
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For each value of one of these n, we get the value of 
x[n]. For clarity, we're often going to color in these 
circles at the top, but that's really just a matter of 
taste. We're either going to label the signal on the y- 
axis or, more typically, in the title of the graph. 

A financial series signal.Daily temperatures over 
the course of a year. The discrete-time plot of a 


speech signal. 


Examples of Discrete-Time Signal Plots 


Here are some examples of signals. The first is a 
financial time series. It's the daily closing share 
price of Google for a five-month period: 


x{[n] 
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You can see here that it's a discrete time signal, 
where each of these signal points corresponds to one 
single share price at the end of a day. There are 
some fluctuations in the price, and if you were a 
financial trader or if you were an economist, you 
would be very interested in the information that this 
daily closing share price signal conveys. Another 
example is a temperature signal, the temperature at 
Houston Intercontinental Airport every day at noon 
for 365 days that comprise the year 2013 (in 
degrees Celsius): 

o. x{n] 
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Again, we can see that there are fluctuations in this 


signal, and if you were a meteorologist or a 
climatologist, you'd be very interested in the 
information that this signal conveys. Finally, here's 
an audio signal that is speech from an actor 


speaking a part in Shakespeare's play, Hamlet: 
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Discrete-time signals are undefined between the 
integer index values and should be plotted 
accordingly. The second plot interpolates between 
the discrete-time integer index values, which is 
inappropriate for a discrete-time plot. 


Plotting Discrete-Time Signals Correctly 


We need to remember that with a discrete-time 
signal, the independent variable is integer valued. 
This means that when you plot a signal in a 
program like MATLAB, you must use a discrete-time 
plotting function (like the stem function) that 
respects the fact that the signal is only defined at 
discrete, integer time points, rather than a function 


(like plot) which interpolates between points: 
x(n] 
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Plotting Complex-valued Signals 


Up to this point, we've been talking about real- 
valued signals. They comprise a single plot of 
n versus x[n]. But what about complex-valued 
signals? 


Recall that a complex number has a real component 
and an imaginary component. There are two 
equivalent ways of expressing a given complex 
number. For some a€C, we can express a in two 
different ways: 


* Cartesian/rectangular form: a= Re(a)+j Im(a) 
* Polar form: a= |alejza, 


where j = —1 (in engineering contexts the variable 

j is used to represent this value because i represents 
electrical current). Just as a complex number can be 
expressed in two different ways, so can a complex- 
valued signal: 


* Cartesian/rectangular 
form: x[n] =Re(x[n]) +j Im(x[n]) 


* Polar form: x[n] = |x[n]|ejzx[n] 


What this means is that, if we're plotting a complex- 
valued signal, we actually need two plots. As we 
have seen, there are two different ways we can plot 
the same complex-valued signal: 


* Cartesian/rectangular form:Plots of the real 
and imaginary parts of a complex-valued signal 
x[n]. 
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* Polar form:Plots of the magnitude and phase of 


a complex-valued signal x[n]. 
|a[n]| 
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Signal Properties 


Signal Classification 


Signals can be broadly classified as discrete-time or 
continuous-time, depending on whether the 
independent variable is integer-valued or real- 
valued. Signals may also be either real-valued or 
complex-valued. We will now consider some of the 
other ways we can classify signals. 


Signal Length: Finite/Infinite 


This classification is just as it sounds. An infinite- 
length discrete-time signal takes values for all time 
indices: all integer values n on the number line from 
— co all the way up to ~. A finite-length signal is 
defined only for a certain range of n, from some N1 
to N2. The signal is not defined outside of that 
range. 

A periodic discrete-time signal. Note how it repeats 
every 8 time units. 


Signal Periodicity 


As the name suggests, periodic signals are those 
that repeat themselves. Mathematically, this means 


that there exists some integer value N for which x[n 
+N] =x[n] , for all values of n. So if we define a 
fundamental period of this particular signal of 
length, like N=8, then we will see the same signal 
values shifted by 8 time indices, by 16, —8,—16, 


etc. Below is an example of a periodic signal: 
x(n 


So periodic signals repeat, and clearly periodic 
signals are going to be, therefore, infinite in length. 
It's also important to remember that to be periodic 
in discrete-time, the period N must be an integer. If 
there is no such integer-valued N for which x[n 
+N] =x[n] (for all values of n), then we classify the 
signal as being aperiodic. 

An infinite-length signal x[n] (only part of it 
shown) has a portion extracted via windowing to 
create a finite-length signal y[n]. An infinite-length 
signal x[n] (only part of it shown) has a portion 
extracted via windowing to create a finite-length 
signal y[n]. A finite-length signal x[n] is periodized 
to create an infinite-length signal y[n] (only a 
portion of it is shown). 


Converting Between Infinite and Finite 
Length 


In different applications, the need will arise to 
convert a signal from infinite-length to finite-length, 
and vice versa. There are many ways this operation 
can be accomplished, but we will consider the most 
common. 


The most straightforward way to create a finite- 
length signal from an infinite-length one is through 
the process of windowing. A windowing operation 
extracts a contiguous portion of an infinite-length 
signal, that portion becoming the new finite-length 
signal. Sometimes a window will also scale the 
smaller portion in a particular way. Below is a 
mathematical expression of windowing (without any 
scaling): 


yt[n]={x[n]N1 <n < N 2 undefined else 


Below is a signal x[n] (assume it is infinite-length, 
with only a part of it shown), with a portion of it 
extracted to create y[n]: 
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There are two ways a signal can be converted from 
finite-length to infinite-length. The first is referred 
to as zero-padding. It is easy to take a finite-length 
signal and then make a larger finite-length signal 
out of it: just extend the time axis. We then have to 
decide what values to put in the new time locations, 
and simply putting 0 at all the new locations is a 
common approach. Here is how it looks, 
mathematically, to create a longer signal y[n] from 
a shorter signal x[n] defined only on N1 <n<N2: 


y[n]={ONOs<n<NIx[n]INI<n<N2 
ON2<n<N3 
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Here, obviously NO<N1<N2<N3, and if we extend 
NO and N3 to negative and positive infinity, 
respectively, then y[n] will end up being infinite- 
length. 


The other way to make an infinite-length signal 
from a finite-length one is to periodize it, which 


means replicating a finite-length signal over and 
over to create an infinite-length periodic version. 
Mathematically, that means defining the new 
infinite-length periodic signal like this: 


y[nJ]=x%m=—-eex[n-mN],NEZ= 
-+x[n+2N]+x[n+N]+x[n]+xI[n 
—-N]+x[n-2N]+°:: 


Graphically, we can see that this amounts to 
repeating the signal over and over, before and after 
the original portion: 

x(n] 
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y(n] with period N = 8 


Periodization and Modular Arithmetic 


It turns out that, as we consider periodization and 
periodic signals, the notion of modular arithmetic 
will be helpful. In modular arithmetic, integers do 
not lie on a line stretching from negative infinity to 
infinity, but rather on a circle of a defined size N. In 


modulo-8, for example, the numbers are 8 spaces on 
a "wheel." Our convention will be for the numbers 
to traverse from 0 to N—1, counter-clockwise: 

O 


Consider a finite-length signal of size N. We can 
align the time-dependent values of the signal on the 
modulo circle: 


When we travel around the clock once, from time 
index 0 to 7, we express the finite-length signal. But 
if we keep traveling, in one direction or the other, 
then that amounts to periodizing the signal. Using 
our modulo notation, we can periodize a finite- 
length signal x[n] to be an infinite-length periodic 
signal y[n] like this:y[n] =x[(n)N]. The modulo 
operation (n)N returns the integer remainder after n 
is divided by N. For example, (13)8 is equal to 5, 


because 8 goes into 13 once, with a remainder of 5. 
For periodization, this would result in 
y[13] =x[(13)8] =x[5]. 


Finite/Periodic Signals Relationship 


We have seen that we can take an N-length finite- 
length signal and periodize it to make an infinite- 
length periodic signal with a period of N. By the 
same token, we can also work in reverse and extract 
one period worth of signal from any periodic signal 
to create a finite-length signal. What this means is 
we can consider periodic signals and finite- 
length signals to be essentially equivalent: we 
can consider just one period of a periodic signal (the 
rest of the signal is redundant, by definition), or 
periodize a finite-length signal. They are two ways 
of looking at the same thing, a phenomenon we will 
often see in our study of signals and systems, and 
we will choose the perspective that best suits our 
needs for particular applications. 

A signal x[n] shifted according to the expression 
x[n—m] where m is positive, and negative. A 
periodic y[n] signal shifted one value in time to the 
right, y[n—1]. 


Shifting Infinite-length Signals 


Given some signal x[n], it will often be necessary 


for us to consider that signal shifted in time. We 
denote such a shift mathematically with an 
expression like x[n —m], where m is some integer. If 
m is greater than zero, then x[n—m] will be just 
like x[n], except it will be shifted to the right by m 
time units. If m is less than zero, it will be shifted to 
the left. Here is who that might look for a couple 
values of m: 


x(n] 
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This type of shifting works the same with periodic 


signals: 
y(n] with period N = 8 
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y(n — 1] with period N = 8 
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Of course, for periodic signals, certain shifts actually 
do not have any effect on them. If a signal repeats 
with a period of N, then shifting that signal by any 
integer multiple of N will yield the original signal. 
Take a look at the signal above, which was shifted 
to the right by a time unit of 1. If we keep on 
shifting it until it is shifted 8 time units the result 
will be identical to the original signal. 

Circularly shifting a signal by m amounts to turning 
it counter-clockwise m steps. Here is the original 
signal x[n] and its shifted version x[(n —1)8]. 


Shifting Finite-Length Signals 


Can finite-length signals be shifted, as well? There 
does not seem to be any reason why not. Suppose 
we have some finite-length signal x[n] (of length N) 
and we define another finite-length signal v[n] to be 
v[n] =x[n-—1]. So we have v[1] =x[0] and 
v[2]=x[1] and so on until v[N —1] =x[N-— 2]. But 
what about v[0], what do we put there? And how 
about x[N—1], where is that supposed to go? We do 
not want to invent information to put in v[0], nor 
lose the information of x[N—1]. An elegant solution 
is to periodize x[n], and then consider the relation 


v[n] =x[n—1]. In this case, we now have a value 
for v[0O]: v[O] =x[—1]. Since x[n] is periodic with 
period N, it also happens that x[—1]=x[N-—1], so 
we do not lose that information. 


This kind of operation, for finite-length signals, is 
called a circular shift, and we can express it 
mathematically with the help of our modular 
arithmetic operator. Circularly shifting a finite- 
length signal x[n] by m time units is expressed as 
x[(n—m)N]. It can also be visualized by turning 
x[n] about the circle on which it resides: 


«(5 
We can time reverse a finite-length signal x[n] with 
the mathematical expression x[(—n)8]. Note here 
how the signal values are traversing clockwise, i.e., 


backwards in time. 


Time Reversing Finite-Length Signals 


For infinite length signals, the transformation of 
reversing the time axis x[ —n] is obvious: just flip 
the signal about n=0. But things are not quite so 
obvious for finite-length signals; if a signal is 
defined for, say, n between 0 and N, then what gets 
flipped across the n=0 from the negative side? 
Once again, it turns out the modular arithmetic 
operator can be called in for help. We reverse the 
time axis, modulo N: x[(—n)N]. Below is an image 
of a finite-length (N =8) signal, time-reversed: 

x(6] 
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x(2| 
A causal signal is 0 for all n<0, whereas an acausal 
signal is 0 for all n=0. 


Signal Causality 


A signal x[n] is causal if x[n] =0 for all n<O, it is 
anti-causal if x[n] =O for all n=O, and it is acausal 
if is neither causal nor anti-causal. 


For an even signal, x[—n] =x[n]. For an odd signal, 
x[—n] = —x[n]. 


Even and Odd Signals 


A signal x[n]is defined as even if x[—n] =x[n], and 
odd if x[ —n] = —x[n]. 
x{n] 
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A digital signal with q=2 bits, so D=22=4 levels. 


Digital Signals 


Digital signals are a special sub-class of discrete- 
time signals. While the independent time variable 
for discrete-time signals is integer-valued, the 
dependent variable (i.e., the value the signal takes 
at any given time) can take on any value. However, 
for digital signals, both of these variables are 
discrete-valued. Rather than take any value on a 
continuum, discrete signals take only a limited 
number of values, or levels. Typically, the number 
of levels is expressed as D= 2q, and each possible 
value of x[n] is represented as a digital code with q 
bits. 

Q(x{n]) 
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Key Discrete-time Test Signals 


In our study of discrete-time signals and signal 
processing, there are five very important signals that 
we will use both to illustrate signal processing 
concepts, and also to probe or test signal processing 
systems: the delta function, the unit step 
function, the unit pulse function, the real 
exponential function, sinusoidal functions, and 
complex exponential functions. We will initially 
consider the first four; sinusoids and complex 
exponentials are particularly important and will be 
treated separately. Each of these signals will be 
introduced as infinite-length signals, but they all 
have straightforward finite-length equivalents. 

The discrete-time delta function. A time-shifted 
discrete-time delta function 6[n—m] , wherem=9. 
Using a time-shifted delta function to isolate a 
sample of the signal x[n]. 


The Discrete-time Delta Function 


The delta function is probably the simplest 
nontrivial signal. It is represented mathematically 
with (no surprise) the Greek letter delta: 5[n]. It 
takes the value 0 for all time points, except at the 
time point n=0 where it peaks to the value 1: 


d5[n] = {1n = 0 0else 
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In a variety of important settings, we will often see 
the delta function shifted by a particular time value. 
The delta function 6[n — mlis 0, except for a peak of 


1 at time m: 
d[n — 9] 
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One of the reasons the shifted delta function is so 
useful is that we can use it to select, or sample, a 
value of another signal at some defined time value. 
Suppose we have some signal x[n], and we would 
like to isolate that signal's value at time m. What we 
can do is multiply that signal by a shifted delta 
signal. We can say y[n] =x[n]6[n—m], but since 
that y[n] will be zero for all n except at n=Mm, it is 
equivalent to express it as y[n] =x[m]6[n—m], 
where now x[m] is no longer a function, but a 
constant. The following figure shows how this 
operation isolates a particular time sample of x[n]: 
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y(n] = x(9]o[n — 9] 
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The unit step function. A shifted step function u[n 


—m] with m=5. A shifted step function can be used 
to zero out all values of a signal before a certain 


time index. 


The Unit Step Function 


The unit step function can be thought of as turning 
on a switch. Usually identified as u[n], it is O for all 


n<0O, and then at n=0 it "switches on" and is 1 for 
all n=0: u[n] ={0n<0O1n=0: 


u{n| 
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As with the delta function, it will also be useful for 


us to shift the step function: 
u{n — 5] 
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And, as you might have guessed, we can use a 
shifted step function in a similar way to the delta 
function by multiplying it with another signal. 
Whereas the delta function selected a single value of 
a certain signal (zeroing out the rest), the step 
function isolates a portion of a signal after a given 
time. Below, a step function is used to zero out all 


the values of x[n] for n<5, keeping the rest: 
x(n] 
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Supposing a signal x[n] were not causal, setting m 
to zero and performing the operation 

x[n]u[n] would zero out all values of x[n] before 
n=O, thereby making the result causal. 

The unit pulse function p[n], here with N1= —5 
and N2=3. 


The Unit Pulse Function 


The unit pulse function p[n] is very similar to the 
unit step function in that it "switches on" from 0 to 
1 at a certain time, but then it also "switches off" at 
a later time. We will say it "switches on" at time N1, 
and "off" at time N2: 


pi[n]={On<NI11N1esnSN20nNFN2 
p{n| 
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Of course, rather than use the above piece-wise 
notation, it is also possible to express the pulse as 
the difference of two step functions: p[n] =u[n 
—N1]-—ul[n—(N2+1)]. 

For a>1, the real exponential function increases 
with time. Here a=1.1. For 0<a<1, the real 
exponential function decreases with time (or we 
could say it increases exponentially as the time 
index decreases). Here a=.9. 


The Real Exponential Function 


Finally, we have the real exponential function, 
which takes a real number a and raises it to the 
power of n, where n is the time index: r[n] =an, 
aER. So at n=0, r[n] =a0, at n=1 it equals a, is a2 
at n=2, and so on. As the name suggests, the signal 
will exponentially increase or decrease over time, 


depending on the value of a. 
r[n] 


Real and Complex Sinusoidal Signals 


Discrete-time real and complex-valued sinusoidal 
signals are an incredibly important signal class in 
the study of discrete-time signals and systems. 
Sinusoidal waves show up in all sorts of science and 
engineering applications, but they are particularly 
relevant for signal processing because (as we will 
see later) they are the foundation of Fourier 
analysis. 

(a) A plot of cos(On). At every point in time, this 
signal takes the value cos(0)=1. (b) A plot of 
sin(On). At every point in time, this signal takes the 
value sin(O) =0. (c) A plot of sin(w4n+ 276). (d) A 
plot of cos(stn). Note how when w=n, as in this 
example, the signal is oscillating as rapidly as 
possible, between —1 and 1 at every single time 
instance. This phenomenon is the opposite of when 
w=0, for which the signal does not oscillate at all. 
So in some sense we can see that 0 is the lowest 
possible frequency w, and x is the highest. Four 
cosine waves with the same frequency, but different 
phases: cos(st6n — 0), cos(st6n — 14), cos(at6n 
— 72) =sin(s6n), and cos(z6n — 2st) = cos(a6n). Note 
how a phase shift of 12 shifts the cosine to be a sine 
wave, and a phase shift of 2x shifts it all the way 
over to where it was before the phase shift. 


Real Valued Sinusoids 


There are two real-valued discrete-time sinusoidal 
wave signals: the sine wave signal and the cosine 
wave signal. They are represented mathematically 
as sin(wn+ ) and cos(wn+ Q). Let's take a look at 
those in more detail. First, as we have seen with 
other discrete-time signals, n is the independent 
variable time index, and it runs from negative 
infinity to infinity. The variable w is known as the 
frequency of the sinusoidal signal, and we will see 
how changing the value of w impacts the rate of the 
signal's oscillation. The variable © is the phase of 
the signal, and changing it will shift the signal left 
or right along the time axis. Finally, the terms sin or 
cos return the corresponding trigonometric value 
of wn+@ for each value of the time index n. Here 
are a few examples of real sinusoidal waves: 
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We saw in the figure above how the frequency w 
influences the rate of the wave's oscillation. The 
other variable in the signal, the phase @, can shift 
the wave backwards and forwards along the time 
axis, without affecting the frequency. Below are 
plots of a cosine wave which all have the same 
frequency, but with a variety of phase shifts (..e., 
different values of ©): 
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A three dimensional visualization of a complex 


sinusoid. Note that this image is continuous-valued, 
whereas a discrete valued version would actually 
appear as points along the line. A complex sinusoid 
plotted according to its magnitude and phase. Note 
that the magnitude of a single complex sinusoid is 
trivial, as |ej(@n+)| =1. A complex sinusoid 
plotted according to its real and imaginary parts. 
These are a cosine, and sine, respectively, which 
follows from Euler's Formula. 


Complex Valued Sinusoids 


So we have reviewed the real sine waves sin and 
cos, and perhaps seeing them in proximity brought 
to mind a very special relationship called Euler's 
Formula: ej8 = cos(9) + jsin(8) (vou may remember 
this from math class with an i instead, but recall 
engineers use that letter for current, and we call the 
imaginary number j). That formula works for any 
particular value of 8, so of course it applies when 
we consider wn-+ 0, as above, which gives us a 
complex valued sinusoid: ej(wn + 6) = cos(an 

+ ©)+jsin(@n+ 0). Let's look at some plots of 
complex sinusoids. Unlike two-dimensional real 
sinusoids (which have a one-dimensional 
independent time variable n and a take a one- 
dimensional value at each time value), complex 
sinusoids are three dimensional: they have the time 
dimension, a real dimension, and an imaginary 
dimension. So they can be visualized as a three 
dimensional helix in space: 


If you were to look at this helix from directly above, 
you would see only the real portion of the helix, and 
it would appear to be a cosine wave. If you looked 
at it from the side, you would see the imaginary 
aspect of it, as a sine wave. The frequency variable 
@ controls how quickly the helix rotates across time 
n, and also the direction: positive values cause it to 
rotate in the counterclockwise manner shown, and 
negative values would result in it rotating 
clockwise. 


While it is illuminating to visualize complex 
simusiods in three dimensions, in practice it is 
actually most common to view them in two, 
separately plotting either the real and the imaginary 
parts with respect to time, or the magnitude and 
phase across time: 
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We'll wrap up our introduction of sinusoids by 

briefly considering the concept of negative-valued 
frequencies. It is easiest to see the difference a 
negative frequency makes, compared to a positive 
frequency of the same magnitude, by expressing it 

all mathematically: ej(/—@)n = e 

—jon = cos(—wn)+jsin(—@n) = cos(on) —jsin(wn). 
So negating the frequency of a complex sinusoid has 
no effect on the real part of the signal (the cosine), 
but it flips the sign of the imaginary part (the sine). 


This operation (preserving the real part, but 
changing the sign of the imaginary part) is also 
known as taking the complex conjugate of the 
signal. So negating the frequency of a complex 
sinusoid is the same thing as taking the complex 
conjugate of it:ej(—o)n = e—jwn = (ejmn)*. 


Why use imaginary numbers? 


Now perhaps you are wondering the point of using 
imaginary numbers. After all, aren't all real world 
signals, well, real-valued? They are indeed, but we 
can simply consider them as the real-part of a 
complex-valued signal. And why go to that trouble? 
There are many good reasons, but here is one to 
start with: exponential functions are much easier to 
work with than trigonometric functions. You can 
easily simplify eaeb into a single term, but you very 
likely would be turning to a table to simplify 
sin(a)cos(b), wouldn't you? 


The Peculiarity of Discrete-Time Sinusoids 


Compared to their continuous-time counterparts 
(those that take a continuous-valued independent 
time variable t), discrete-time sinusoidal signals 
have two unique characteristics. It is possible for 
them to alias, and they are not always periodic. 
Here x1[n] =cos(x6n) and x2[n] =cos(1326n) . 
While they have different frequencies, they are 
identical signals, for cos(13st6n) = cos((16 + 2st)n). 
Low frequencies are those w close to 0 or 2x rad, as 
in the first signal, cos(st10n). High frequencies are 
those w close to m or —2z rad, as in the second 
signal, cos(9710n). 


Aliasing of Discrete-time Sinusoids 


One might think that if two different discrete-time 
sinusoids had different frequencies, then they would 
be different signals. Such is always the case with 
continuous-time sinusoids, but not always for the 
discrete-time version. Consider two discrete-time 
sinuoids x1[n] and x2[n] with different frequencies, 
@ and w+ 2m: 


xl[n]=ej(mn+ 0) 
x2[n] =ej((moa+2nx)n+90) 


We can then simplify the expression of x2[n], using 


the fact that ej21=1 to arrive at this surprising 
conclusion: x2[n] = ej((@ + 2m)n+ d =ej(wn + 2an 
+ 0) =ej(wn + Od)ej2an = ej(mn + O)(1)n= ej(wn 

+ 0)=x1[n] 


So x1[n] and x2[n] have different frequencies, yet 
they are identical! You can see this plotted out with 
w=s6 below: 
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This phenomenon is called aliasing. It happens 
when frequencies are offset by any integer multiple 
of 2% (you can use w and w+2zk in the example 
above to see for yourself). 


So only frequencies along a continuous interval of 
length 2 on the real number are distinct from each 
other. For this reason, when we deal with discrete- 
time frequencies we consider only those along the 
interval O<w <2 or —1<O <Q, as any other 
frequency aliases back to an identical signal with a 
frequency in that range. Within these ranges, 
frequencies close to 0 (or 21, depending on the 
range used) are low frequencies--their sinusoids do 


not oscillate very quickly-- and frequencies close to 
m (or —st, depending on the range) are high 
frequencies: 
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x1[n] =cos(27316n) is periodic, with N= 16. 
x2[n] =cos(1.16n). The frequency of 1.16 is not a 
fraction of 2st, so this sinusoidal signal is not 
periodic. 


Periodicity of Discrete-time Sinusoids 


Recall that a signal x[n] is defined to be periodic if 
there exists some integer N for which x[n 
+N]=x[n] , Vn. Continuous-time sinusoids are 
always periodic, but again, such is only sometimes 
the case for discrete-time sinusoids. Suppose we 
have a complex sinusoid with a frequency w = 2xkN, 
where k and N are integers. This just means that w 
is some fraction of 2x. It turns out that this signal is 
periodic, with period N: 


x[n] = ej(2akNn + 0)x[n + N] =ej(2akN(n 
+N)+0)=ej(2akNn + 2nkNN + 0) = ej(2akNn 


+ O)ej(2nkNN) = ej(2akNn 
+ O)ej(2akNN) = ej(2akNn + ©)(ej(20k)) =x[n]. 


Here is a plot of a sinusoid with frequency 27316. 
You will note that it has a period of N=16: 
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Now, these fractions of 21 are special values of 
we call harmonic frequencies, for sinusoids with 
such frequencies are periodic. 


In contrast, consider sinusoids whose frequencies 
are not fractions of x: 


x[n] =ej(wn+o)x[n+N] =ej(w(n 

+N)+0)=ej(mon+oN+0)=ej(o@n 

+ O)ej(woN)=ej(won+ dO)ej(oN) = xl 
n],unless oN = 2nk~0=22kN 


So we see that discrete-time sinusoids are periodic 
if, and only if, their frequencies are fractions of 2s. 
Consider the example of a non-periodic sinusoid 
below. It definitely oscillates, and at first it appears 
to be periodic, but look carefully and you will see 
that it actually is not, unlike the one from the figure 
above. Not only is this cosine not periodic, but aside 
from at respective negated indices (e.g., —1 and 1, 
—2 and 2), it will never even take the same value 


more than once! 
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Take note of sinusoids whose frequencies are of the 
form w= 2zkN, for they will play a starring role in 
the Fourier analysis of periodic and finite-length 
discrete-time signals. 


Complex Exponentials 


From Complex Sinusoids to Complex 
Exponentials 


Recall the form of a discrete-time complex sinusoid: 
x[n] = ej(wn+ 0). As we have already seen, that 
signal itself is complex-valued, i.e., it has both a real 
and an imaginary part. But look closely at just the 
exponent, and you will see that the exponent itself is 
purely imaginary. 


Suppose we let the exponent be complex-valued, say 
of the form a+jb, where a and Db are real. Then we 
have e(a+jb)n =eanejbn = (ea)nejbn. The result is a 
complex sinusoid multiplied by a real exponential 
signal (whose base is ea). 

The real and imaginary parts of a complex 
exponential xn for which |z| <1. The real and 
imaginary parts of a complex exponential xn for 
which |z|>1. 


Complex Exponentials, Defined 


We do not typically represent complex exponentials 
in the way derived above, but rather express them 
in the form x[n] =zn, where z is a complex number. 
Being a complex number, it lies on the complex 


plane with a magnitude of |z|, and an angle of zz 
we define as w. So then, if we would like to express 
x[n] =zn as a combination of a real exponential and 
a complex sinusoid, as above, we have: x[n] =zn=| 
z|nejown. Below are some plots of complex 


exponentials for different ae of z. 
Re(z"), |z| <1 


-15 -10 -4 0 5 10 15 


on 
Im(z”), |z| <1 


Tt, ‘iat ot tte eettteesecsseoees 


-15 -10 “5 0 ; 0 15 


Re(z Nae |>1 
eeeessttonseetttegsgetttag 
‘15 -10 -5 0 5 10 5 

rn 

Im(z"), |z| > 1 


tteevenentneegeetttes get te 


-15 -10 “5 0 5 10 15 
r 


We see that when the magnitude |z| is greater than 
1, the signal oscillates and exponentially grows with 
time, and if the magnitude is less than 1, it decays 
over time. If the magnitude is exactly equal to 1, it 


does not grow or decay, but only oscillates. In fact, 
if the magnitude is 1, the complex exponential is, by 
definition, simply a complex sinusoid: |1| 
nejwn=ejon. Therefore you can see that complex 
sinusoids are a subset of the more general complex 
exponential signals. They correspond to values of z 
on the complex plane that lie on the unit circle, |z| 
=1. 


Signals are Vectors 


One mathematical way of understanding signals is 
to see them as functions. A signal x[n] carries some 
kind of information, having a value x[n] at every 
given point of its independent time variable n. 
Another, and complementary, way of understanding 
signals is to consider them as vectors within vector 
spaces. By doing this we will be able to apply 
various tools of linear algebra to help us better 
understand signals and the systems that modify 
them. 


Vector Spaces 


A vector space V is a collection of vectors such that 
if x,y€V and a is a scalar then axEV and x+yEV 
(along with some other properties we would expect, 
such as commutativity, additive and multiplicative 
identities, etc.). In words, this means that if two 
vectors are elements of a vector space, any 
combination or scaled version of them is also in the 
space. When we consider the scaling factors, we 
mean a that are either real or complex numbers. 


There are many different kinds of vector spaces, but 
the two in which we are especially interested are RN 
and CN. RN is the set of all vectors of dimension N, 
in which every entry of the vector is a real number, 


and CN is exactly the same, except the entries are 
complex valued. 


Starting Small - A Two Dimensional 
Vector Space 


You are already familiar with a prominent example 
of a vector space, the two-dimensional real 
coordinate space R2. Every ordered combination of 
two real numbers is a vector in this space, and can 
be visualized as a point or arrow in the two- 
dimensional Cartesian plane. Suppose x and y are 
each vectors in R2. x=[x[0]x[1]] , y=[y[0ly[1]] 
where x[0],x[1],y[0],y[1]€R. The indices we use to 
refer to elements of the vector start their numbering 
at 0. This is the common convention in signal 
processing and many programming languages, like 
C, but note that vector indices in MATLAB start with 
1. Scaled versions of these vectors are still within 
the space: 


ax=a[x[O0]x[1]]=[ax[O]ax[1]]€ 
R2 


So is the sum of two vectors: x 


+y = [xLO]x{1]] + LylOly{1]] = [x[0] + y[O]x[1] + yl1}JER: 


When a vector x in RN is scaled by a&R, the result 
is still in RN. Here is an example of an N- 
dimensional signal/vector x, which is then scaled by 
a=3. The sum of two vectors in RN is another 


vector in RN. 


The Vector Space RN 


Let's now move from two dimensional vectors to 
those with N dimensions, each taking a real valued 
number: 


x= [x[0]x[1]:x[N—1]], 


where x[n]€R So we have a mathematical entity x, 
with N ordered real values associated with it. Stated 
this way, we see that x is a signal, simply expressed 
in a vector form x as opposed to the signal notation 
form x[n], but both forms refer to the exact same 
thing. 


Just as with in 2-dimensions, we can perform the 
same operations on the N-dimensional signal/vector 


x(n] 


) 5 10 15 20 25 1] 


When each of the real-valued elements of x is scaled 
by a real value a, the results are still real-valued, of 
course, so the resulting scaled vector is obviously 
still within RN. As RN is a vector space, that should 
come as no surprise, and neither should the fact that 
the sum of two vectors in RN is itself another vector 
in the vector space RN: 
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The Vector Space The Vector Space CN 


Ordered sequences of N complex numbers form the 
vector space CN, just like their real counterparts do 
in RN: x= [x[0]x[1]:x[N—1]], 


where x[n]€C . So each element of a vector in CN is 
a complex number, which can be represented either 
in terms of its real and imaginary parts, or in terms 
of its magnitude and phase. The vector as a whole 
can also be represented in those ways, either in 
rectangular form 


x=[Re xLO0]+j Im x[O]Re x[1]+j 
lor xi Re oe LN 1] ey ko N= 
i= Rees PaO Voki LIN ap el 
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or in polar form 


x = [|x[0]|ejzx[0] |x[1]|ejzx[1]:|x[N—1]|ejzx[N 
—1]]. 


Linear Combinations of Vectors 


In the study of signals and signal processing, there is 
a particular mathematical operation that will show 
up quite a few times, that of a linear combination. 
Given a collection of vectors in a vector space, say 
M vectors x0,x1,...xM—1€CN, and M scalars a0,a1, 
...,aM —1€C, then the linear combination of these 
vectors is: 


y=a0x0+a1x1+-:-+aM—-1xM-1=2m=0M 
—lamxm. 


The result is itself also a vector in the vector space 
(indeed, this is one of the requirements of what it 
means for something to be a vector space). 


Linear Combination Example: A Mixing 
Board 


A linear combination is a scaled sum of different 
vectors in a vector space. A real world example of a 
linear combination is the mixing board used in 
music recording studios. The board takes in a 
variety of different inputs and combines them-- 
amplifying some, reducing the level of others--to 
create a single output of music. Mathematically, we 
could say x0 is drums, x1 is bass, x2 is guitar, ..., 
x22 is a saxophone, and x23 a singer, and then the 


mixing board creates the linear combination 
y=a0x0+ a1x1+-:-+aM—-1xM-1=2m=0M 

— lamxm. Changing the different am's would result 
in a different kind of sound that either emphasizes 
or de-emphasizes certain instruments. For example, 
the producer may be particularly interested in 
increasing the volume of the cowbell. 


Linear Combinations as Matrix Multiplication 


It is possible to express a linear combination 
without explicitly using sums, but rather as the 
product of a matrix and a vector. To do this, all of 
the vectors to be scaled and summed in the linear 
combination must first be arranged into a matrix: 

X = [xO|x1]---|xM — 1] = [xO[0]x1[0]---xM 
—1[0]x0[1]x1[1]---xM —1[1]:::xO[N —1]x1[N—-1]::: 
xM—-1[N-1]] 


Then the scaling factors for each vector are stacked 
into a single vector: 


a=[a0al:aM-1] 


Once those are in place, the linear combination is 
ultimately expressed as a simple matrix-vector 
multiplication:y = a0x0 + a1x1+-:-+aM—1xM 
—1=%m=0M- 1lamxm = [x0|x1|---|xM—1] 
[a0al:aM — 1] = Xa. 


The Strength of a Vector 


In the study of signal processing, the "strength" of 
signals (which, remember, can be expressed as 
vectors) is often of interest. It is not immediately 
clear what standard could be used for judging the 
strength of a vector. A similar problem comes up 
when determining the physical strength of a person: 
does it have to do with the amount of weight that 
can be lifted, or the number of times a 

certain weight can be lifted, or perhaps how far the 
weight can be thrown? So it is with vectors; 
consider the two below: 


Which signal is "stronger"? 
a| nj 
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Which is stronger than the other? Perhaps you 
might have intuition just by looking at them, but 
undoubtedly it will be better if we can introduce 
some kind of rigorous definition of strength. 


Norms: The Strength of Vectors 


There indeed are such ways to measure the strength 
of vectors, and they are called norms. Just as there 
are different ways to measure human strength, there 
are also different ways to measure the strength of 
the vector. The first one we will consider is the 2- 
norm. 


The Euclidean length, or 2-norm, of a vector xECN 
is defined as: 


Ix]2=2n=ON-1/x[n]]|2 


This norm, like all norms, is a function that takes a 
vector as its argument and then produces a real 
number that is always greater than or equal to 0. 
The 2-norm is a very common norm, so often it will 
be seen without even having a subscript, i.e. as |[xll. 


Let's calculate a 2-norm. Suppose: 

pao el ee 

Then ||x||2, the 2-norm of x, is: 

IIxI2 =n = ON — 1|x[n]|2=124+ 224+ 32=14. 

It happens that the 2-norm is actually just one kind 


of a more general type of norm, called p-norms. For 
any real number p=1, the p-norm is: 


Ixllp=(Zn=ON-—1|x[n]|p)1/p 


So the 2-norm is a p-norm with p=2. Another norm 
of interest to us is when p=1, the 1-norm: 


Ixl1 =2n=ON-1]|x[n] | 


A final norm that we will consider is called the ~- 
norm. It is essentially what happens as the p in the 
p-norm increases to ~, but has this easy to calculate 
expression: ||x||ec = maxn|x[n]]|. So if you have a 
vector x and take the absolute value of all the 
entries, the ~ norm of x is simply the largest of 
those values. 


The Meaning of Different Norms 


There is a reason for having different norms; they 
each are helpful in their own way for measuring 
particular things about a signal. Some of these 
measurements have clear physical quantities 
associated with them. For example, consider the 
electrical signal travelling through a loudspeaker, 
and suppose x corresponds to the voltage of this 
signal. The 2-norm of this signal is significant, 
because it corresponds to the electrical energy 
going through the coil (although this is actually a 
bit of a simplification). The ~-norm of the signal is 
important because it represents the maximum 
voltage value going through the circuit. We might 


want to pay attention to each of those values, 
though for different reasons. If the 2-norm is too 
large, then the excessive amount of energy might 
melt the circuit (or lead to a large electric bill!). If 
the ~-norm is too large, the extreme voltage at a 
given time might cause the speaker cone to move 
too far, breaking it. 


Or consider another example of the utility of 
different norms. Suppose you are designing an 
autopilot system for a car. Among other concerns, 
you clearly will want to make sure it does not 
unnecessarily deviate from its lane. Call the 
deviation at a given time the error signal d. Would 
you seek to minimize ||d||2 or |Idlloo? If you minimize 
Id|[2, then the car will, overall, keep to its lane very 
closely, but quick and large deviations would be 
permitted: if over the course of a long drive the car 
stayed in its lane perfectly, except for a fraction of a 
second it accidentally swerved into another lane, the 
Id|[2 would be small, but that could potentially still 
be a disaster! Instead, you would want to minimize 
the ||d|loo value, so the car never moves more than a 
maximum distance from, say, the center of the lane. 


Normalizing a Vector 


Suppose you run a recording studio and you are 
putting different recorded tracks onto a single 
album. You run in to the following problem: the 


artist was standing closer to the microphone on 
some tracks, and farther away on others. The 
loudness varies wildly across the different tracks, so 
it would be frustrating for customers to listen to the 
album, having to change the volume with every 
track. What would you do? 


If you have studied signal processing, then you will 
know that you can consider each track as a signal, 
and that each track will have a different norm. In 
order to make the loudness the same, you would 
like to modify each track so that it has the same 
"strength", or norm (technically, unless all the tracks 
are the same length, you would like them to have 
the same norm-squared per time unit; and to get even 
more technical, you would also need to make other 
modifications based upon the frequencies of the 
input). Thankfully, it is very straightforward to 
modify a signal so as to change its norm. There are 
two steps. 


First, you will modify the vector so that it is a unit 
norm, i.e. anorm of 1. This is called normalizing 
the vector. To do that, you simply scale the vector 
by its norm. Suppose your vector is x; the 
normalized version of x is then x|[x||2. For example, 
earlier we saw that the vector: 


aS [123 


had a norm of 14. The normalized version of x then 


is x divided by that value: 
x= [114214314] 


That new vector has a 2-norm of 1, as desired, and 
is oriented in the exact same direction as the 
original vector. If you wanted the norm of the 
vector to be some other value, say the number M, 
you would multiply the normalized vector by M. 


Inner Products and Orthogonality 


In discrete-time signal processing, understanding 
signals as vectors within a vector space allows us to 
use tools of analysis and linear algebra to examine 
signal properties. One of the properties we may 
want to consider is the similarity of (or difference 
between) two vectors. A mathematical tool that 
provides insight into this is the inner product. 


Transposing Vectors 


One of the ways to express the inner product of two 
vectors is through matrix multiplication. To explain 
that, we must first introduce the concept of 
transposing vectors. In order to multiply two 
matrices (a vector is simply a matrix in which one of 
the dimensions is 1), the column dimension of the 
first must match the row dimension of the second. 
To make those match for two vectors with the same 
dimensions, we must transpose one of them. To 
take the transpose of a matrix, simply turn the rows 
into columns: the first row will become the first 
column in the transposed matrix; the second row, 
the second column, and so on. Here is how that 
looks for a vector. An N-row, single column vector 
transposed becomes a 1 row, N-column vector: 


[x[O]x[1]:x[N — 1]]T=[xL0]x[1]---x[N—-1]]. 


Now, when it comes to complex-valued vectors, we 
can take a transpose in the same way, but for the 
purposes of finding an inner product we actually 
need to take the conjugate, or Hermitian, transpose, 
which involves taking the transpose and then the 
complex conjugate of each entry: 


Lx Oe ee EN ee ta Ee od Ed 
eox[N—1]*] 


Of course, for real valued vectors, the regular 
transpose and Hermitian transpose are identical. 


The Inner Product 


The inner product of two complex (or real) valued 
vectors is defined as: 


(x,y) =y Bx = Tn] 0N=— ix bal yin 


So the inner product operation takes two vectors as 
inputs and produces a single number. It turns out 
that the number it produces is related to the angle 0 
between the two vectors: 


cos(@x,y)=Re (x,y)Ilxll2lly|]2 
This formula works for complex and real vectors, 


although taking the real part of the inner product is 
redundant for real-valued ones. 


For two- (or three-) dimensional vectors, this angle 
is exactly what you would expect it to be. Let 


x=[12], y=[32] 


We have |[x|[22=12+ 22=5, |lylI]22=32+22=13 
and (x,y) =(1)(3) + (2)(2) =7. The angle between 
them is arccos(7(5)(13)) =0.519 rad = 29.7°. If you 
plot the vectors out in the Cartesian plane, you will 
indeed see an angle between them of about 30 
degrees. 


For higher dimensional signals the result of the 
inner product--how it relates to the angle between 
signals--may not seem as intuitive, but the 
information it provides is still just as useful, and of 
course it is computed in the same was as with 


shorter vectors. Consider signals below: 
& in| 
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The inner product of these two signals, computed 
according to the formula above, is 
(x,y) = yTx = 5.995, which corresponds to an 


angle of 0x,y = 64.9°. 

The inner product of these two signals is 0. Despite 
overlapping nonzero entries in the time domain, the 
inner product of these two particular signals is also 
0. 


Inner product: Limiting Cases 


Recall that the inner product is defined as a sum 
(xn = ON — 1x[n]y[n]*), which can also be expressed 
as a vector product (yHx). Let's look at a couple of 
interesting values that sum could take. 


The dot product of two signals could be rather large. 
If the signals are identical, it is simply the norm of 
the signal, squared: 


(x,x) = Zn = ON — 1x[n]x[n]* = Xn=ON — 1|x[n]| 
2=|Ix|l22. 


On the other hand, it is also possible for the dot 
product sum to be 0. Consider the two signals 
below: 


le 
O.5 + 
ou 


0 5 10 15 20) 
n 


The inner product of those two signals is obviously 
0 because each point-wise product is also 0. But it is 
possible, of course, for products in the sum to be 
nonzero and still have the total add up to 0: 


0 5 10 


20 25 


int. 


) 5 10 15 20 25 


rn 


Whenever the inner product of two signals is 0, it is 
defined that those signals are orthogonal. 


Orthogonality of Harmonic Sinusoids 


Recall the special class of discrete-time finite length 
signals called harmonic sinusoids: 

sk[n] =ej2aNkn, n,k,NEZ, O<n<N-1, O<k<N 
ab 


It is a very interesting property that any two of 


these sinusoids having different frequencies (i.e., 
k#]) are orthogonal: 


(sk|sl) =Xn=ON-1sk[n]sl*[n] =2 
n=ON-lej2naNkn(ej2aNln)*=2n 
=ON-lej2mNkne—-j2xaNln=2Xn=0 
N-lej2aN(k-—1)n lett r=k-1€Z,r 
=O=%n=ON-1lej2nNrn=2%n=O0N - 
lan with a=ej2xNr,andrecall Xn = 
ON =Lan = 1 =—:aN1l-= a= 1-> ¢).2076N.N 
l—ej2nrN = 0 V 


If two of these sinusoids have the same frequency, 
then their inner product is simply N: 


(sk,sk) =|Iskl22=z2n=ON-1|sk[n] 
|2=xyn=ON-1lej2nkNn|2=2Xn=0 
N-11= Nv 


So the inner product of two harmonic sinusoids is 
zero if their frequencies are different, and N if they 
are the same. In order to make the latter number 1, 
instead of N, they are sometimes normalized: 


s~k[n]=1Nej2xNkn, nik N-S 7 
Osn=N-1, OsksN-1 


Matrix Multiplication and Inner Products 


Let's take look at the formula for the matrix 


multiplication y = Xa. For notation, we will 
represent the value on the nth row and mth column 
of X as xm[n]. The matrix multiplication looks like 
this: 


Va Saley (al) = Pee lar er ad 
M-1l1[n]l:::]l[aQ0al:aM-1] 


And the value for the multiplication is: 

y[n] =xm =0M — lamxm[n]. Hopefully that formula 
looks familiar! What it is showing is that the matrix 
multiplication y = Xa is simply the inner product of 
a with each row of X, when X and a are real. If they 
are complex-valued, then the complex conjugate of 
one of them would have to be performed before the 
matrix multiplication for each value of y to be the 
inner product of the matrix row with a. 


The Cauchy Schwarz Inequality 


Above we saw that the inner product of two vectors 
can be as small as O, in which case the vectors are 
orthogonal, or it can be large, such as when the two 
vectors are identical (and the inner product is 
simply the norm of the vector, squared). It turns out 
that there is a very significant inequality that 
explains these two cases. It is called the Cauchy 
Schwarz Inequality, which states that for two 
vectors x and y, |(x,y)| <Ilxllllyll. Now the magnitude 
of the inner product is always greater than or equal 


to 0 (being 0 if the vectors are orthogonal), so we 
can expand the inequality thus: 0 < |(x,y)| <llallllyll. If 
we divide the equation by ||xllllyll, then we have 0<| 
(x,y) |IIxllllyll <1. This explains why we can define 
cos(Ox,y) = Re (x,y) |Ixll2Ilyll2, for the cosine function 
also has a range of 0 to 1. 


Now there are many different proofs of the 
inequality, and it is something of a mathematical 
pastime to appreciate their various constructions. 
But for signal processing purposes, we are more 
interested in the utility of the inequality. What it 
basically says is that--when the norms of two vectors 
are taken into consideration--their inner product 
ranges in value from 0 to 1. Because of this, we can 
see that the inner product introduces some kind of 
comparison between two different vectors. It is at its 
smallest when they are, in a sense, very different 
from each other, or orthogonal. It is at its peak 
when the vectors are simply scalar multiples of each 
other, or in other words, are very alike. 


It turns out there are many applications in which we 
would like to determine how similar one signal is to 
another. How does a digital communication system 
decide whether the signal corresponding to a "0" 
was transmitted or the signal corresponding to a 
"1"? How does a radar or sonar system find targets 
in the signal it receives after transmitting a pulse? 
How does a computer vision system find faces in 
images? For each of these questions, the similarity/ 


dissimilarity bounds established by the Cauchy 
Schwarz inequality help us to determine the answer. 


Norms and Inner Products of Infinite Length Vectors 


As with finite length signals (which can be 
understood as vectors in RN or CN), we can take the 
norms and inner products of infinite length signals, 
as well. For the most part, the concepts will apply 
just the same. 


The Norms of Infinite Length Vectors 


As with finite-length vectors, the p-norm for infinite- 
length vectors is a sum: ||xllp = (2n = — © o|x[n]| 
p)1/p, p=1. Unlike the p-norms of finite-length 
signals, the summation for infinite-length vectors is 
not necessarily finite. As a consequence, the p- 
norms of infinite length signals may not be 
bounded. Consider a finite-length signal x[n] =1 for 
0<n<N-—1. The 2-norm for this signal would be N. 
However, for an infinite length signal of constant 
value 1, all the p-norms will be unbounded! 


For some infinite-length vectors, the p-norm may 
exist only for certain values of p. Consider, for 
example, the signal: 


ent HO n =O Lan "1 


0 : 10 lf 0 25 


The 2-norm of this signal exists: 


[|xI[22 = Xn = — oo co |x[n]|2=Xn=1-|1n| 
2=%n=1°1n2=726 = 1.64. 


But the 1-norm is unbounded: ||x||1 = Zn = — © | 
x[n]|=Xn=1leln=, 


For a finite-length vector, the ~-norm is simply the 
maximum value of the magnitudes of all its 
elements. For infinite-length vectors, it is instead the 
supremum, or the least upper bound of the set of 
all the elements of the vector. The distinction 
between that value and the maximum value of the 
set is outside the bounds of this course. For all the 
cases we will consider, the --norm will simply be 
the largest of the magnitudes of the elements of the 
vector. For example, for the signal x[n] above, 
lIx[n]llec = 1. 


Inner Product of Infinite Length Signals 


The inner product of infinite length signals is 
defined the same way as for finite-length signals, 
except that the sum is infinite: 


(x,y) =Zn= — e cx[n]y[n]*. The relationship to the 
angle between the signals is also the same: 
cosOx,y =Re (x, y)I|Ixll2lly||2. 


Linear Combinations of Infinite Length 
Vectors 


When it comes to infinite-length vectors, we will 
also be interested in linear combinations of them. 
However, in contrast with finite-length vectors, we 
will be interested in the linear combination of an 
infinite number of the vectors: y= Xm = 

— co coqmxm. As with finite-length vectors, infinite- 
length vectors can be stacked into a matrix (an 
infinitely large one) and be multiplied by an infinite 
length vector of scalers to perform linear 
combination: 


XS | eed | KO ed 


a=[:a-—la0dal:] 

y==im=—-erwoeamxme= [i:ixk-1[- 
1]xO[-—1]x1l[-1]-::: x-—1[0]x0O[0] 
Lu. On| ee x-1[1]xO[1]x1[1]-":::l[ 


:a-—-ladOal:]=Xa 


Discrete-time Systems 


Discrete-time signals are mathematical entities; in 
particular, they are functions with an independent 
time variable and a dependent variable that 
typically represents some kind of real-world 
quantity of interest. But as interesting as a signal 
may be on its own, engineers usually want to do 
something to it. This kind of action is what discrete- 
time systems are all about. A discrete-time system 
is a mathematical transformation that maps a 
discrete-time input signal (usually designated x) to a 
discrete-time output signal (usually designated y). In 
other words, it takes an input signal and modifies it 
to produce an output signal: 

System H takes takes a discrete time signal x as an 
input and produces an output y. 


Ab H Y 


There is no end to the possibilities of what a system 
could do. Systems might be trivially dull (e.g., 
produce an output of 0 regardless of the input) or 
incredibly complex (e.g., isolate a single voice 
speaking in a crowd). Here are a few examples of 
systems: 


* A speech recognition system converts acoustic 
waves of speech into text 


¢ A radar system transforms the received radar 
pulse to estimate the position and velocity of 
targets 

¢ A functional magnetic resonance imaging 
(fMRI) system transforms measurements of 
electron spin into voxel-by-voxel estimates of 
brain activity 

* A 30 day moving average smooths out the day- 
to-day variability in a stock price 


Signal Length and Systems 


Recall that discrete-time signals can be broadly 
divided into two classes based upon their length: 
they are either infinite-length or finite-length (and 
recall also that periodic signals, though infinite in 
length, can be viewed as finite-length signals when 
we take a single period into account). Likewise, 
discrete-time systems are also finite- or infinite- 
length, depending on the kind of input signals they 
take. Finite-length systems take in a finite-length 
input and produce a finite-length output (of the 
same length), with infinite-length systems doing the 
same for infinite-length signals. 


Examples of Discrete-time Systems 


So a system takes an input signal x and produces an 
output signal y. How does this look, 
mathematically? Below are several examples of 
systems and their mathematical expression: 


Identity: y[n] =x[n] 

Scaling: y[n] = 2x[n] 

Offset: y[n] =x[n]+ 2 

Square signal: y[n] =(x[n])2 

Shift: y[n] =x[n+m],mEZ 

Decimate: y[n] =x[2n] 

Square-time: y[n] =x[n2] 

Moving average (combines shift, sum, and 
scale): y[n] = 12(x[n] + x[n—1]) 
Recursive average: y[n] =x[n] +ay[n—1] 


So systems take input signals and produce output 
signals. We have seen some examples of systems, 
and have also introduced a broad categorization of 
systems as either operating on finite- or infinite 
length signals. 


Linear Systems 


Discrete-time systems are mathematical 
transformations that take input signals and map 
them to output signals. For a given input x, a 
discrete-time system will produce a new signal y: 


4 H Y 


It turns out that it is very important to know or 
determine if a system has certain characteristics, 
and among these is the classification of linearity. A 
system is linear if it has two special 

properties: scaling (sometimes called homogeneity) 
and additivity. The first of these is satisfied if, for 
any arbitrary input x, scaling the input by any 
complex valued constant value a will result in the 
output being scaled by the same value. 
Mathematically, this scaling rule is represented as 
H{ax} = aH{x}, and as a system diagram it looks 
like this: 

A system that satisfies the scaling property. 


ee 


A system has the additivity property if, for any two 
arbitrary inputs, the output of the sum of them is 
the same as the sum of their individual outputs: 
H{xl +x2}=H{x1}+H{x2}. 
A system that satisfies the additivity property. 

x1 Yr x2 Y2 


A system has both the additivity and scaling 
properties (and thus, by definition, is linear) if for 
arbitrary inputs x1 and x2 and arbitrary constants 
al and a2, H{alxl1 +a2x2}=alH{x1} + a2H{x2}. 
If a system lacks either of the scaling or additivity 
property, then it is said to be nonlinear. 


Determining Linearity 


Consider again the definition of linearity; to be 
linear, a system must preserve the scaling and 
additivity properties for any arbitrary input. 
Therefore, determining linearity amounts to 
completing a mathematical proof that assumes 
arbitrary inputs and scaling factors, and then shows 
the necessary result. Suppose H{x[n]} = 3x[n]. Here 
is how the additivity proof would look: 


Let x1 and x2 be arbitrary inputs to system H and 
let a1 and a2 be arbitrary constants. 


H{x1[n]} = 3x1 [n]H{x2[n]} 

= 3x2[n]H{alx1[n] + a2x2[n]} 

= 3(a1x1[n] +a2x2[n]) = 3a1x1[n] + 3a2x2[n] = al1(3x1 [1 
— Linear. 


Now, to show a system is nonlinear requires a 
different kind of proof. Rather than having to prove 
both of the properties hold for any arbitrary 
input(s), only a single example needs to be provided 
for which either of the properties fail. For example, 
consider the system H{x[n]} =x[n] +1. We can 
show it is nonlinear by considering just the single 
following case: 

Let x[n] =OH{x[n]}=x[n]+1 1 
= 1 But 42s Pn} S=2 x [a] eS 
1=1+2+2H{x[n]}— Nonlinear 


=0+ 
2-0 + 


Good students of signals and systems must 
become adept at determining the linearity (or 
nonlinearity) of systems. Practice on the system 
examples below; which of them are linear? 


¢ Identity: y[n] =x[n] 

* Scaling: y[n] = 2x[n] 

* Offset: y[n] =x[n] +2 

¢ Square signal: y[n] = (x[n])2 

¢ Shift: y[n] =x[n+m]mEZ 

* Decimate: y[n] =x[2n] 

¢ Square time: y[n] =x[n2] 

* Moving average (combines shift, sum, 


scale): y[n] =12(x[n] +x[n—1]) 
* Recursive average: y[n] =x[n]+ay[n-1] 


Identity: Linear 

Scaling: Linear 

Offset: Nonlinear 

Square signal: Nonlinear 
Shift: Linear 

Decimate: Linear 

Square time: Linear 
Moving average: Linear 
Recursive average: Linear 


OOO O0O0000 


Linear Systems and Matrix Multiplication 


One of the interesting characteristics of any linear 
system is that its input/output relationship can be 
expressed as a matrix multiplication (note that this 
is distinct from--though ultimately related to--the 
notion of using matrix multiplication to express one 
signal as a linear combination of others). In fact, 
this relationship is actually an identity: any linear 
system can be expressed as a matrix multiplication, 
and matrix multiplications are linear systems. Below 
is how to represent any linear system 
mathematically, either in matrix multiplication 
notation, y = Hx or, in summation 


notation, y[n] = Xm[H]n,mx[m] = Xmhn,mx[m] 
(where hn,m = [H]n,m represents the row-n, 
column-m entry of the matrix H). 


This matrix multiplication can be understood in two 
ways. First, the multiplication means that each 
value in the vector y is the inner product of the 
corresponding row of H with the vector x. Or, 
equivalently, the vector y can be seen as a weighted 
sum of the columns of H, with the values in the 
vector x being the weights of the corresponding 
columns. Below is a picture of matrix multiplication, 
with different colors representing different values. 
Try to comprehend the multiplication with both 
perspectives. 


Time-Invariant Systems 


Recall that a discrete-time system is a mathematical 
entity that takes an input signal (usually denoted x) 
and produces an output signal (usually denoted y). 
In the study of signal processing, systems that have 
certain characteristics are of particular interest. 
Among these characteristics is time-invariance. 


There is a very basic intuition behind the notion of 
time-invariance. In many cases, we would like a 
system to behave a certain way, no matter when an 
input may be given. To use a practical analogy, 
people expect their toasters to operate the same way 
on Tuesdays as they do on Mondays. 


A system is time-invariant if a time delay for an 
input results in the same time delay on the output. 
Expressing this idea of system time-invariance 
mathematically is straightforward. Consider a 
system H; let x[n] be some arbitrary input, and call 
the system's output for that input y[n]. H is time- 
invariant if, for some arbitrary integer value q, 
H[x[n—q]]=y[n—q]: 


x(n yin] 
z[n — q| y[n — q] 


Example: Moving Average System 


Consider, for example, the moving average system 
y[n] =12(x[n] +x[n—1]). Is this system time 
invariant? To find out, let's delay the input by some 
value q, and see what the output is. We'll call the 
delayed input signal x’[n] (=x[n—q]) and the new 
output y’: 


y [n] = 12(%’[n] + x’[n—-1]) =12(«%[n—-q]+x[(n 
=1) = ql): 


Now the important question, is y’[n] =y[n—q]? 
That will be easy to determine, simply replace every 
n with n—q in that original y[n] equation and see if 
it is the same as y’[n]:y[n—q] =12(x[n—q]+x[n-—q 
~1))=12(«[n—q] +xl[@—-1)-q)=y'Inlv 


Example: Decimation 


Now let's consider a decimation system, 

y[n] =x[2n]. Is this system time-invariant? First, 
create a new signal that is a delayed version of x[n]: 
x’[n] =x[n—q]. Next, find the output corresponding 
to this new signal: y’[n] =x’[2n]. Now express this 
output in terms of the original input: y’[n] =x 

‘[2n] =x[2n—q]. Finally, check to see if this is the 


same as y[n—gq]. As y[n—q] =x[2(n—q)] is not 
equivalent to y’[n] (since 2n—1+42(n—q) for all n 
and q), the system is not time-invariant. 


Now practice on the system examples below; 
which of them are time-invariant? 


¢ Identity: y[n] =x[n] 

* Scaling: y[n] =2x[n] 

* Offset: y[n] =x[n] +2 

¢ Square signal: y[n] = (x[n])2 

¢ Shift: y[n] =x[n+m]mEZ 

¢ Square time: y[n] =x[n2] 

* Recursive average: y[n] =x[n] +ay[n-1] 


¢ Identity: time-invariant 

¢ Scaling: time-invariant 

¢ Offset: time-invariant 

¢ Square signal: time-invariant 

¢ Shift: time-invariant 

¢ Square time: not time-invariant 

¢ Recursive average: y[n] =x[n] +ay[n 
— 1]time-invariant 


Time-Invariance for Finite-Length 
Systems 


All of the discussion above has considered systems 
which take infinite-length signals as inputs. A 
system which takes finite-length signals as inputs 
and produces them as outputs can also exhibit the 
characteristic of time-invariance, albeit with a 
slightly different definition. Finite-length systems 
are time-invariant if any circular shift on the input 
results in a corresponding circular shift on the 


output: 
a yl 
x(n aw ul(n— a) 


Linear Time-Invariant Systems 


Among the many characteristics and classifications 
of discrete-time systems, two of particular 
importance are linearity and time-invariance. If a 
system happens to exhibit both of these qualities, 
then it is referred to as being an LTI system (linear 
time-invariant). These systems are very significant 
in the study of signal processing, for reasons that 
will be clear when the concept of the system 
impulse response is considered. Using the tests for 
both linearity and time-invariance, determine 
whether or not the following systems are LTI. 


Consider the systems below. Which are LTI? 


* Identity: y[n] =x[n] 

* Scaling: y[n] =2x[n] 

* Offset: y[n] =x[n] +2 

¢ Square signal: y[n] = (x[n])2 

¢ Shift: y[n] =x[n+m]mEZ 

* Decimate: y[n] =x[2n] 

¢ Square time: y[n] =x[n2] 

* Moving average (combines shift, sum, 
scale): y[n] =12(x[n] +x[n—1]) 

* Recursive average: y[n] =x[n] +ay[n-1] 


¢ Identity: LTI 
* Scaling: LTI 


* Offset: not LTI (time-invariant, but not 
linear) 

¢ Square signal: not LTI (time-invariant, 
but not linear) 

¢ Shift: LTI 

¢ Decimate: not LTI (linear, but not time- 
invariant) 

¢ Square time: not LTI (time-invariant, but 
not linear) 

* Moving average (combines shift, sum, 
scale): LTI 

* Recursive average: LTI 


A pictorial representation of a Toeplitz matrix. 


Matrix Structure of LTI Systems (infinite- 
length) 


We recall that linear systems can be expressed 
through a matrix multiplication. Suppose that a 
system is linear; any output y can be expressed as 
the multiplication of an infinite-dimensional matrix 
H with the input x: 


y = Hxy[n] = ~Xm[H]n, mx([m] 


Now if a linear system is also time-invariant, it turns 
out the matrix H will have an interesting structure. 
To see this, we will first express the matrix 
multiplication in a summation notation, where 


hn,m = [H]n,m represents the row-n, column-m 
entry of the matrix H: 


y(n] = H{x[n]} = XEm=-—-eowhn, 
mx[m],-—-° <n< o 


Supposing the system is time-invariant, we have: 
H{x[n—q]} = Xm= — © ohn,mx[m-—gq] = y[n-—q] 


If we apply a simple change of variables (n’=n—q 
and m’=m-—dq), we then have: 


H{x[n’]} = Xm’= — © ohn’+q,m’+qx[m’] = y[n 
1. 


If we compare this final equation with the original 
one, we see that for an LTI system, hnm = hn+q,m 
+ qVqEZ. 


So for LTI systems, the matrix H that defines the 
system's input/output relationship has a special 
structure: 


hn,m = hn+qg,m+qVqEZ 


(where hn,m is the value at the nth row and mth 
column of the matrix H). Such matrices are called 
Toeplitz matrices. They have identical values along 
their diagonals: 


H = [iii«h—-1,-1h—-1,0h-1,1-- 
-~h0O,—-1h0,0h0,1°"-* h1,—1h1,0h1 


pleeiii] = [iii-ho0,0h—1,0h- 2,0 
ieee h1,0h0, Tae eee me a 
One) 


Because of that property, all of the information in 
the matrix is contained in a single row, or column. 
We'll call the Oth row h[n]: h[n] =hn,0: 


H = [iii-+-h0,0h-—1,0h-2,0----h1, 
Oh0,Oh-—1,0-+---h2,0h1,0h0,0-°33:] 
= [ili-h[OJh[—1]h[-2]---h[1] 
h[OJh[—1]--h[2]h[1]h[O] 2:2] 


Below is an example of a pictorial representation of 
a Toeplitz matrix. Note how the diagonals are all 
the same color: 


“4 2 0 2 4 
™m 
A pictorial representation of a circulant matrix. 


Matrix Structure of LTI Systems (finite- 
length) 


For systems operating on finite-length signals, we 
can apply similar analysis. Linear systems can be 
expressed as a matrix multiplication y = Hx, here 
shown using summation notation: 

yin] = H{x[n]} = Xm=ON-—1hn,mx[m],0<n<N 
—1. If we enforce time-invariance on such a 
system, we then have: 


H{x[(n—q)N]} = =2m=ON-— 1hn,mx[(m 
—Q@N] = y[m—q)N]. 


A change of variables (n’=n—q) then yields: 


H{x[(n’)N]} = XYm’=—-qM-1-gh 
(n’+q)N,(m’+q)Nx[(m’)N] = yl 
n’)N] 


We see, then, that for LTI finite-length systems, 
hnm = h(n+q)N,(m+q)NVq€Z. This is similar to 
the Toeplitz structure for infinite-length systems, 
but note the circular shifts, how the rows "wrap" 
around the edges of the matrix: 


H = [h0,0h0,1h0,2:°*h0O,N—-1h1,0 
his thio hi.N = Lh?..0ho. Tho. 
ho, N=} 2bN TORN =—1, LAN <1; 
2--hN—-1,N—1] = [h0,0hN-1,0h 
N—2,0°-h1,0h1,0h0,0hN-1,0°h2 
,0h2,0h1,0h0,0°h3,0:!::hN—1,0h 
N-2,0hN-3,0°h0,0] 


This kind of matrix is called a circulant matrix. As 
with the Toeplitz matrices, all of the information in 
a circulant matrix is contained in a single row or 
column, the other rows/columns simply being 
(circularly) shifted versions of the original. So we 
can call the N-length signal h[n] the Oth column of 
the matrix, and express the whole matrix in terms of 
it: 


H = [h0O,OhN-1,0hN-2,0°h1,0h 
1,0h0,0OhN—-1,0°-h2,0h2,0h1,0h0 
,0h3,0::::HN-1,0OhN-2,0hN-3, 
O--h0O,0] = [h[O]h[N-1]h[N-2] 
~hHE1T]h[1]h[OJh[N—-1]-h[2]h[2] 
HiLIWlOl“WiSls thin =1TIhIN =21 
h[N-—3]-h[0]] 


Below is a pictorial representation of a circulant 
matrix. Like Toeplitz matrices, the values along 
diagonals are equivalent, but in addition to having 
this property, the values also "wrap around" the 
matrix boundaries: 


The Impulse Response of Discrete-Time Systems 
The system H takes an input x and produces an 
output y. The system matrix of an LTI system has a 
Toeplitz structure. Note how when the H matrix is 
multiplied by it, the delta function "selects" the Oth 
column. A pictorial representation of the impulse 
response and system matrix of a finite-length LTI 
system.A recursive average system 

(y[n] =x[n] + ay[n—1]) is an example of an IIR 
system, as the impulse response has an infinite 
number of nonzero values. The impulse response of 
the moving average system above 

(y[n] = 12(x[n] +x[n—1])) has only 2 nonzero 
values, and thus is an example of an FIR system. 


The Impulse Response of Infinite-Length 
LTI Systems 


Systems are mathematical transformations that take 
input signals and map them to output signals: 


r H y 


Recall that linear systems can be expressed as 
matrix multiplications, the output y being the 
product of the infinitely long matrix H and the input 
vector x: y=Hx. Recall also that if a linear system is 
also time-invariant, the matrix H has a special 


"Toeplitz" structure: each column/row is simply a 
shifted version of any other column/row. 


The Oth column of the H matrix for an LTI system 
has a special interpretation. If we multiply the H 
matrix by a delta function, then the result will 
simply be that Oth column: 


H —_ rs 


We call the result of this multiplication, i.e. the 
multiplication of the H matrix by delta (impulse) 
function, the impulse response of the system. It is 
usually denoted h[n], or if referring to it as a vector, 
as h. The impulse response h[n] is the output of a 
system when the input is 6[n]. Using the matrix 
notation from above, h is the Oth column of the 
system matrix H. 


Because the impulse response is simply a column 
(specifically, the Oth column) of the system matrix 
H, then if H is an LTI system, the impulse response 
completely characterizes the system. It tells us 
absolutely everything we need to know about the 
system! Why is this? Because we can construct the 
entire system matrix H by simply shifting the 
impulse response accordingly. The value of the 
matrix at row n and column m is h[n—m]. 


As a result, we can also express the system input- 
output matrix multiplication in terms of the impulse 
response. If we have that [H]n,m is the value of the 
H matrix at row n and column m, then the matrix 
multiplication y[n] =Hx=%m = — © o[H]n,mx[m]. 
From above we have the matrix values expressed in 
terms of the impulse response; subbing this in, we 
see the output is a sum involving the impulse 
response and the input: y[n] =Hx=2Xm= — © oh[n 
—m|]x[m]. 


Impulse Response of the Scaling System 

Let's find the impulse response of a simple scaling 
system: y[n] =H{x[n]} = 2x[n]. To find the impulse 
response h[n], simply calculate the output when 
the input is a delta function: 

h[n] = H{d[n]}) = 28[n]. 


h[n] = 2 4[n] 


If we would like to know what the system matrix H 
is, we use the formula [H]n,sm=h[n-—m] =26[n 
—m]. Below is a pictorial representation of that 
matrix: 

The system matrix of a simple scaling system. 


Impulse Response of the Shift System 


Now let's try finding the impulse response of a 
simple time-shifting system: y[n] =H{x[n]} =x[n 
— 2]. Once again, find the impulse response by 
supplying a delta function as the input: 

h[n] =H{S[n]} =8[n—2]. 


h{n] = d[n — 2] 
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The system matrix H can be found, as before, as 
[HJn,m=h[n-—m]=s[n—m-—2]. 
The system matrix of a shift system. 


Impulse Response of a Moving Average System 
Consider now a moving average system: 

yin] = H{x[n]} = 12(x[n]+x[n-—1)]). For this 
system, we find that the impulse response is: 

h[n] = H{$[n]} = 12(6[n] +3[n—1)). 


h{n| = 4 (5[n] ; d[n 1}) 
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The system matrix is found to be [H]n,m=h[n 
—m]=12(6[n—m]+6[n—m-—1]). 
The system matrix of a moving average system. 


Impulse Response of the Recursive Average 
System 

recursive moving average system has the 
following input-output 
relationship:y[n] = H{x[n]} = x[n]+ay[n—1]. 
To find the impulse response of this system, we 
first let x[n] =6[n] and find the 
output: h[n] =H{x[n]}+ah[n—1]. A bit more 
work is required if we would like a non-recursively 
defined impulse response. Assuming that the 
system is initially at rest (it has no output in the 
absence of an input), then we can find the output 
by computing several values of h[n]: 
h[ — 1] =Oh[0] =8[0] +ah[—1] =1+a0-0=1h[1] = 41] +<¢ 
Having discerned a pattern, we see that 
h[n] =anu[n] (this may be more rigorously proven 
by induction, if desired): 


h[n] = a" u[n], a=0.8 
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rn 
pictorial representation of the recursive average 
system. 


The Impulse Response of Finite-Length LTI 
Systems 


For a finite-length LTI system, just as with an 
infinite-length system, each row/column of its 
system matrix is a shifted version of another row/ 
column. Thus, the system matrix for a finite-length 
system is also of a Toeplitz structure. However, for a 
finite-length system, not only are the rows/columns 
shifted, they are circularly shifted; the values in the 
rows/columns "wrap around" when shifted across 


the boundary of the matrix. 


As a result, the entire system matrix can be 
expressed in terms of the Oth column of the matrix, 
which (as with infinite-length systems) is given the 
mathematical expression h[n] and is named the 
impulse response. For a length-N LTI system, the 
formula to determine the value of the matrix at a 
particular row n and column m is [H]n,m = h[(n 
—m)N]. Note how this is nearly identical to the 
infinite length formula, except for the modulo-N 
operator. 


As with infinite-length systems, being able to 
express the entire system matrix in terms of the 
impulse response means that the system's input- 
output matrix multiplication relationship can be 
expressed as a sum involving the input and impulse 
response: 


y(n] =Hx=XYm=ON-1h[(n-m)N]x 
[m]. 


Impulse Response Length of LTI Systems 


For LTI systems, it is often of interest whether its 
impulse response is of infinite length (i.e., infinite 
nonzero support) or finite length. Such a distinction 
is typically made only for infinite-length systems (as 
the impulse response of finite-length systems is 
obviously always finite-length!). If a system has an 
impulse response of infinite length, it is called an 
infinite impulse response (IIR) system; if its 
impulse response is of finite length, it is called a 


finite impulse response (FIR) system. 
h[n] =a" u[n], a=0.8 
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h{n] = 4 (d[n| + d[n — 1}) 


Convolution of Infinite-Length Discrete-Time Signals 


It is a fundamental operation in signal processing, 
given an LTI system H and input x[n], to find the 
output y[n]. From our study of LTI systems, we have 
seen that there are several ways to find the output 
of the system: 


* One option is to find the output y[n] at each 
time n using the input-output relationship 
(which defines the output y[n] in terms of the 
input and output at other various times). 

¢ Another option is to take the system matrix H 
and the input vector x and perform a matrix 
multiplication y = Hx. 

¢ A final option is to find the output at a given 
time in terms of a summation involving only 
the input and the system's impulse 
response: y[n] = Xm= — © oh[n—m]x[m]. 


That final option has a special name: it is the 
convolution of x[n] and the system impulse 
response h[n]. It also has a special operator symbol, 
*, so that the convolution of x[n] and h[n] is 
expressed as x[n]*h[n]. 


Computing Convolution 


Given an LTI system's impulse response h[n] and an 


input x[n], the process of computing the output y[n] 
via the convolution sum formula (Xm = — © eoh[n 
—m]x[m]) is straightforward, according to the 
following steps: 


* Step 1: decide which of x[n] or h[n] you will 
flip and shift; you have a choice, since 
x[n]*h[n] =h[n]*x[n]. We will suppose you 
choose to flip h[n]. 

* Step 2: Plot x[m]. 

* Step 3: Plot h[—m], the time-reversed version 
of the impulse response. 

* Step 4: To compute y[n] at the time point n, 
plot the time-reversed impulse response after it 
has been shifted to the right (i.e, delayed) by n 
time units: h[ —(m—n)]=h[n-m]. 

* Step 5: y[n] is the inner product between the 
signals x[m] and h[n—m] (technically it is the 
sum of the point-wise products; for complex 
signals, do not complex conjugate the second 
signal, as you would for a true inner product). 

* Step 6: Repeat for all n of interest (potentially 
all n€Z). 

* Step 7: Plot y[n] and perform a reality check to 
make sure your answer seems reasonable. 


Here we will illustrate convolution by convolving 
two signals together, x[n] and h[n], plotted below. 
e will call the result of the convolution y[n], so 


yin] =x[n]*h[n]. 


rn 
Recall that the first step of the convolution process 
is to decide which of the signals we will "flip and 
shift." The final result will be the same in either 
case, the only difference is in working the problem 
out. But here it truly does not matter, because the 
signals are identical. But we will say that we are 
flipping and shifting h[n]. 
The second and third steps are to plot x[m] and 
li Eton) 


Now step 4 and 5 are to shift h[n—m] by varying n 
and then computing the inner product between h[n 
—m] and x[m]. For each shift, that inner product 
will be the value of the output y[n] for that value 
of n. In the figures below, we will see how h[n—m] 
and x[m] look for each n (x[m] of course will not 
change), and will progressively build up y[n] as n 
increases. We'll start with n= — 1: 


x[m| 
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Notice how for this value of n, and for all n less 
than that as well (for which h[n—m] will be even 
farther to the left), x[m] and h[n—m] have no 
overlapping nonzero values. Therefore the inner 
product in these cases is zero. 
So for n<0O, the dot products of x[m] and h[n—m] 
are zero, and hence y[n] is zero for these n. Now 
let's consider n= 0: 
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So now, for n=0, x[m] and h[n—m] overlap just a 
bit (at m=0), and the value of the dot product of 
the two signals is 1. So y[0] = 1. 
The next shift to the right makes n=1: 
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There are two nonzero locations where x[m] and 
h[n—m] line up, m=0 and m=1. The dot product 
of the signals for this n=1 is 1:1+1-:1=2, so 
y[2] =2. 
Continuing to shift h[n—m] to the right, we have 
fOr: 


x{m| 
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For n=2, the two signals are now maximally 
aligned, with a dot product of 1:1+1-:1+1:1=3. So 
y[2] =3. 
t n=3, h[n—m] has moved past the point of 
maximum alignment; as with n=1, the dot product 


between it and x[m] is 2: 


x[m| 
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We're approaching the end of our work on this 

convolution problem. At n=4, h[n—m] has just 

about moved on past the point of overlapping 

monzero values with x[m]: 


y(n] 
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The dot product between the two signals is only 1. 
Finally, for n=5, there is no nonzero overlap 
between x[m] and h[n—m], so the dot product is 
zero. For all n>5 (for which h[n—m] moves even 
farther to the right) this is the case as well: 


7" z{m| 
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We have found the dot product between x[m] and 
h[n—m], and therefore y[n], for all integer values 
of n, so our convolution computation is finished. 
The final step of the convolution is to examine the 
plot of y[n] and do a quick "reality check" to make 
sure the answer is reasonable. At this introductory 


point, you might not know what to expect a 
convolution looks like, so it may not be clear if it is 
reasonable or not. But in time, especially as you 
understand the "flip and shift" nature of 
convolution, you will gain a feel for how 
convolutions should generally look. 


Convolution of Finite-Length Discrete-Time Signals 


Given a discrete-time, finite-length LTI system with 
an impulse response h[n], it is a common operation 
to find the system's output y[n] for some input x[n]. 
There are several ways this can be done, including 
finding the system matrix H and then using it to find 
the output via matrix multiplication: y = Hx. 
Recalling the circulant structure of that matrix, the 
output can be expressed in formula 

form: y[n] = x[{n]@h[n] = Xm=ON-1h[(M 
—m)N]x[m]. This operation is known as circular 
convolution (or finite-length convolution). 


Computing Circular Convolution 


The circular convolution x[n]®h[n] can be 
computed using the following seven step procedure: 


1. Decide which of x[n] or h[n] you will flip and 
shift. Since circular convolution is 
commutative, either choice will result in the 
same output. 

2. Plot the values of x[m] on a "wheel," beginning 
on the right side of it (for m=0), and 
continuing counter-clockwise for a total 
number of places that is equal to the length of 
the signal. 

3. Plot h[{(—m)N] on another wheel. This simply 


amounts to plotting the values in a clockwise 
direction. For a given value of n, y[n] will be 
the dot product of the two "wheels." Your 
original plots of x and h correspond to n=0. 
For n=1, first rotate the values of hina 
counter-clockwise manner, i.e., delay them one 
time unit. Rotating h in this way is what 
creates h[(n—m)N]. Then... 

. Rotate the h wheel accordingly and... 

. Perform the inner product. 

. Repeat this for n ranging from 0 to N—1. 

. Plot your final answer and do a "reality check" 
to see if it seems reasonable. 
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Let's see the circular convolution process in action. 
Consider the two finite-length (N = 8) signals 
below: 

Discrete-time finite-length signals x[m] and h[m], 
N=8. 


x{n| 
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Since N=8, our "wheel" will have have 8 places. 


We will place the initial value on the right side of 
the wheel, but the starting place is arbitrary, so 
long as you are consistent (you can also choose to 
have your signal run clockwise, if you prefer). 

For the purposes of computing circular convolution 
by hand, the values of the signals may be plotted 
on a "wheel" with a number of places equal to the 
length of the signals. The time values progress 
counter-clockwise. 


The second step is to plot x[m] and h[(— m)8]. For 
x, this means starting on the right and then writing 
the values counter-clockwise. For h, this means 
writing the values clockwise: 

Plots of x[m] and h[(n—m)8]. The values of x[m] 
are plotted counter-clockwise. To plot h[(—m)8], 
the values of h[m] are simply plotted clockwise, 
i.e., backwards in time, starting at the same point. 


Having plotted x[m] and h[(—m)8], we can now 
begin to find our output values corresponding to 
each n. We'll start with n=0, which for h[(n—m)8] 
is of course what we have originally plotted: 

Plots of x[m] and h[(0—m)8], the inner products 
of which is y[0]. The plot of x[m] will remain 
constant for each of the inner products. For the 
plot h[(n—m)8], the h[(—m)8] signal is rotated n 


places units counter-clockwise. Here it is plotted 
for n=O. 


We'll move around counter-clockwise, starting at 
zero position, to find the inner product. The value 
of the inner product for n=0 is: 

y[0] =1:0+2-0+ 
—1:0+0-0+0-4+0-3+0-2+0-1=0. 

Now, for each successive value of n, starting with 


n=1 and going up to n=7, we will rotate the h 
plot counter-clockwise one time unit, and then find 
the inner product with the (unchanged) x. 

x[m] and h[((n—m)8] for n=1. 


The value of the inner product for n=1 is: 

y[1] =1:1+2-:0+ 
—1:00+0:0+1:0+0-4+0:3+0:2=1. 

Below are the plots for n= 2. x remains unchanged, 


and h is shifted one time unit counter-clockwise 
from where it was at n=1: 
x[m] and h[(n—m)8] for n=2. 


The value of the inner product for n=2 is: 
viZi=11—20- 
—1:0+00+0:0+0-4+0:3+0-2=2+2=4. 
Now for n=3: 

x[m] and h[(m—m)8] for n=3. 


The value of the inner product for n=3 is: 
y[3] =1:34+2:24+ 
—11+00+0:0+00+00+0-4=3+4-1=6. 
Next, n= 4: 

x[m] and h[((n—m)8] for n=4. 


The value of the inner product for n= 4 is: 
y[4]=1:44+2:34+ 
—1:2+0:1+00+00+00+00=4+6—-2=8. 
For n=5: 

x[m] and h[((n—m)8] for n=5. 


The value of the inner product for n=5 is: 
y([5] =1:0+2-44 
—1:3+0-2+0:1+0-0+0:0+0:0=8—3=5. 
For n=6: 

x[m] and h[(—m)8] for n=6. 


The value of the inner product for n=6 is: 

y[6] =1:0+2:0+ —1-4+0-34+0-2+0:1+00+0-0= 
nd finally, n=7: 

x[m] and h[(n—m)8] for n=7. 


The value of the inner product for n=7 is: 
y[7] =1-:0+2-0+ 

—0+04--0:3--02--01 -00—0. 

Now that we have found the output for each n 
from 0 to 7, we may plot y[n]. 

The final plot of y[n]. 


Properties of Discrete-Time Convolution 


Recall the definition of discrete-time convolution. 
For infinite-length signals we 

have: y[n] = x[n]*h[n] = Xm= — ~ oh[n 
—m]x[m],-— 9° <n<o, 


The formula for finite-length (with length N) signals 
1S: 


y[n] = x[n]@h[n] = Xm=ON-1h[(n 
—m)N]x[m],0<n<N-1. 


As with other mathematical operations, convolution 
has several important properties which we will now 
explore. We will express these properties in terms of 
infinite-length convolution; unless otherwise noted, 
they apply also to finite-time convolution. 

The commutativity property of convolution, as 
displayed in block diagram form. 


Convolution is Commutative 


Like addition and multiplication, the operation of 
convolution is commutative. For two signals x[n] 
and h[n], x[n]*h[n] =h[n]*x[n]. This means that 
when convolution is calculated by hand, either of 
the two signals can be the one chosen to "flip and 
shift." It also means that the output of a system with 
impulse response h[n] and input x[n] is the same as 


a system with impulse response x[n] and input h[n]. 
In block diagram form, these two diagrams are 
equivalent: 


Though the output remains the same, appreciating 
the commutativity of convolution may give added 
insight to the convolution process. 


Proving the commutativity of convolution is 
straightforward; simply apply a change of variable 
to the convolution sum equation: 

y(n] =x[n]*h[n] = Xm = — ~ oh[n 


—m]x[m]Let v=n-—m—~m=n-v=2n-v=-— ©-n 
—v= oh[v]x[n-—v] =Xv= © +nv= — © + nh[v]x[n 
—v] =Yv= — ov= oh[v]x[n—-v] =Xv= 


— ov=ox[n—vjh[v] =h[n]*x[n] 

A "cascade" of LTI systems. A cascade of two LTI 
systems is equivalent to a single LTI system with an 
impulse response that is the convolution of the 
cascaded systems' impulse responses. 


Convolution is Associative: Cascade of 
LTI Systems 


Suppose a signal is processed by one LTI system, 
and the output is fed into another LTI system: 
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Can these two "cascaded" systems be simplified so 
that they are represented by a single system? It 
turns out that they can, because convolution is not 
only commutative, it is also associative: 

(x*h1)*h2 =x*(h1*h2). Therefore if a signal x[n] is 
input into a system with impulse response h1[n], 
and the output of that system is then the input a 
system with impulse response h2[n], the final 
output y[n] is equivalent to that if x[n] had been 
input into a signal system with an impulse response 
of h1*h2: 
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Proving the associativity of convolution is a matter 
of careful rearrangement of the sums: 
(x*h1)*h2[n] =*m= — © oh2[n—m](1= 

— co cooht[m—1]x[]]) =X1]= — © oXM= — w~ oh2[n 
—mJh1l[m—1]x[]] =x] = — © o(Ym= — © oh2[n 
—mJh1l[m—1])x[]] =x]= — ~ © (h1*h2)[n 

—])x[1] =x*(h1*h2)[n]. 

A parallel combination of LTI systems. A parallel 
combination of two LTI systems is equivalent to a 
single system whose impulse response is the sum of 
the parallel systems' impulse responses. 


Convolution is Distributive: Parallel LTI 
Systems 


Suppose now that one input x[n] is fed into two 
different LTI systems, and then the two outputs are 
then summed to form a new output: 


| a | - 
m= 
Just as with the cascaded systems, this parallel 
arrangement of LTI systems can also be expressed as 
a single system, because of the distributive property 


of convolution, that 
x[n]*h1[n] +x[n]*h2[n] =x[n]*(h1[n] + h2[n]): 


a 


The proof of this property is a consequence of the 
distributive property of 

addition:x[n]*h1[n] +x[n]*h2[n] = Xm= — © ohl1[n 
—m]x[m] + *m= — © oh2[n—m]x[m] =xym= 

— co coo(hl[n—m]x[m] +h2[n—m]x[m])=*m= 
—occo(hl[n—m]+h2[n 

—m])x[m] =x[n]*(h1[n] +h2[n]). 


Discrete-Time Infinite-Length and Finite-Length 
Convolution Equivalence 

When a signal x[n] of duration length of 4 is 
convolved with itself, the duration length of the 
result is 4+4—-1=7. 


The Nonzero Support of the Convolution 
of Infinite-Length Signals 


The convolution of two infinite-length signals yields 
another infinite-length signal. However, just because 
a signal is infinite in length does not mean it has an 
infinite number of nonzero values. The delta 
function 5[n] is nonzero only at n=0, even though 
it is infinite-length (i.e., it is defined for all integers 
Nn). 


For reasons that will be apparent shortly, it is worth 
considering the duration interval of infinite length 
signals and their convolutions. Suppose an infinite- 
length signal is nonzero at all points less than n=0 
and greater than n=5. Its duration interval goes 
from n=0 to n=5, and we say that the length of 
this duration is 6. It turns out that if signal a[n] has 
a duration length of Da, and signal b[n] has a 
duration length of Db, than the duration length of 
the convolution a[n]*b[n] is Da+Db—1. 


"Equivalence" of Linear and Circular 
Convolution 


The reason why we would like to consider the 
duration length of two infinite-length signals' 
(linear) convolution is that there are situations in 
which we would like to produce the same result 
through circular convolution (i.e., the convolution 
of finite-length signals). Of course, the outputs of a 
linear convolution and circular convolution are not 
equivalent, strictly speaking; they produce two 
distinct categories of signals: the first one infinite- 
length, the second one finite-length. It is possible, 
however, to take the duration sections of infinite- 
length signals and treat them as finite-length signals 
so that their (circular) convolution is the same as the 
duration section of the infinite (linear) convolution of 
the infinite-length signals. 

The finite-length signal x~[n] here is simply the 


duration section of the infinite-length signal x[n] 
from above. The circular convolution x~[n]@x~ [n] 
is clearly not the same as the duration section of the 
convolution of x[n] with itself. Once again, the 
circular convolution in this case is clearly not the 
same as the nonzero support of the convolution of 
x[n] with itself. 


Zero-padding and Resulting Circular 
Convolutions 


This is all perhaps best explained with illustrations. 
Consider the following infinite-length signal x[n] 


(obviously only a portion of it is shown): 
z{n] 


This convolution of this signal with itself produces 
another infinite-length signal (again, only a portion 
of it is shown): 

x[n] * x[n] 


We will attempt to produce finite-length signal that 
is equal to the duration section of that convolution 


above (the nonzero portion, from n=0 to n=6). To 
do this, we will first try to circularly convolve two 
finite-length signals that are simply the duration 


section of x[n]. We will call this signal x~[n]: 
&[n 


So it is apparent that simply circularly convolving 
the duration sections of x[n] will not do. We 
consider what would happen if we add a 0 to the 


end of x~[n] and try again. 
z\n 


Well, this does seem to be a step in the right 


direction of having x~[n]@x~[n] equal the 
duration section of x[n]*x[n]. Let's add another 0 to 
x~[n]: 

x{n] 


&[n] @z [n] 


We are getting close now, and perhaps you have 
already figured out how many zeros we need to 
ultimately add to achieve the desired circular 
convolution. Let's add another: 


z{n] 


&[n] @ &{n) 


So it seems that three zeros added to make x~ [n] 
did the trick, for the circular convolution 
x~[n]@x~[n] is indeed equal to the duration of 


x[n]*x[n]. Let's see, though, what happens if we add 
yet another zero: 


z[n] 


At this point, it seems that adding extra zeroes to 
x~[n] will also add zeroes to the circular 
convolution; there's no harm with that, but it isn't 
necessary for our goal of the circular convolution 
equaling the duration section of the linear 
convolution. 


Requirements of Circular/Linear 
Convolution Equivalence 


Suppose we have two signals that produce some 
result via linear convolution; if we would like to 
achieve that same result (or at least, the duration 
section of that result) through circular convolution, 
then we will need to cut out the relevant sections of 
the two signals and then zero-pad them before 
circularly convolving them. 


Let's define the amount of zero-padding with a bit 
more rigor. Consider infinite-length discrete-time 
signals a[n] and b[n]. Let anz[n] and bnz[n] be 
finite-length signals which are the duration sections 
of a[n] and b[n], respectively; that is, they are the 
smallest contiguous sections of a[n] and b[n] 
outside of which all values in time before and after 
those sections are zero (it is, of course, allowable 
that some values within the contiguous sections may 
be zero). Call the duration length of these duration 
sections Da and Db, respectively. Now zero-pad the 
signals anz[n] and bnz[n] (i.e., appended zeroes to 
the end of the signals) to bring the new length of 
these signals to be Da+ Db—1, and call the new 
zero-padded signals a~n and b~n. Then we have 
that a~n®b~n is equivalent to the duration section 
of a[n]*b[n]. 


Put more simply, suppose that the duration length 
of a[n]*b[n] is D (which we have seen equals Da 
+Db-—1). We must zero-pad the duration sections 
of a[b] and b[n] to each be a length of D to have 
their circular convolution be equivalent to the 
duration section of the linear convolution. 


Why this All Matters 


There is a very good reason to go into this 
somewhat tedious task dealing with duration 
sections and their lengths and the relationships of 


linear and circular convolution. The reason has two 
parts. First: real-world input/output LTI systems 
implement linear convolution. Second: there are 
tremendous computational complexity advantages 
to performing circular convolution (which we will 
see later). So, if we can express the output of linear 
convolution through circular convolution, it means 
we can perform real-world system computations 
quickly. 


Impulse Response and LTI System Causality 


In addition to linearity and time-invariance, there 
are other significant classifications of discrete-time 
systems. One of these is causality. A system is 
causal if its output, for any n, depends only on 
inputs at or before time n. Causality is an important 
practical qualification on a system: it is not possible 
to implement a real-time system whose current 
output depends on future values! While the 
qualification of causality is also relevant for non-LTI 
systems, it has a special relationship with LTI 
systems. Recall that for LTI systems, the impulse 
response can be used to find the system's output 
given some input (through convolution of the input 
and the impulse response). Likewise, for LTI 
systems, the impulse response can also tell us 
whether or not the system is causal. 

An example of an impulse response for a causal 
system. The impulse response is zero for n<0O. 


The Impulse Response and Causality 


There is a straightforward relationship between an 
LTI system's impulse response and whether or not 
the system is causal: An LTI system is causal if and 
only if its impulse response is O for all n<O (..e., the 
impulse response is a causal signal). This follows 
naturally from the convolution sum. The system's 
output is defined as y[n] = Xm= — © oh[n 


—m]x[m]. Note that if and only if 

h[m] =0, V m<0O, no x[m] for m>n contribute to 
the sum, which is to say that no future values of the 
input contribute to the output at time n. Below is an 
example of the impulse response for a causal 
system. Note that it is O for n<0O. This also 
corresponds to the system matrix being lower 
triangular: 


h[n] = a" u[n], a=0.8 


Impulse Response and LTI System Stability 


It is of practical significance in the design of 
discrete-time systems that they be "well behaved," 
meaning that for any "well behaved" input, the 
system gives a "well behaved" output. Colloquially 
speaking, we do not want an innocuous input to 
result in the system "blowing up." The technical 
term for "well behaved" systems is that they are 
stable. It is essential for many systems that they be 
stable, for the sake of safety and proper operation 
within wider systems. For example, with steering or 
braking or aircraft control systems, it could be 
catastrophic if a small input led to a wildly 
divergent output. 


Consider the recursive average system 

y[n] =x[n] + ay[n—1], with an eminently 
reasonable and contained input of the step function 
u[n]. For values of a less than 1, the system is "well 
behaved," but the output "blows up" for a>1: 
When the input to a recursive average system is a 
step function, the output may or may not be 
contained, depending on the value of the a 
coefficient. For a<1, the system's output is 
contained, whereas for a>1, the system's output 
increases with time exponentially. 


xn] = u[n] 


n 
y(n], with a= 4 


4 2 0 2 ; é 8 
n 
y[n], with a= 3 


Now as with any desirable characteristic or quality, 
there are many ways to define stability. We could 
say a system is stable if its output never exceeds a 
particular value, or perhaps that the output's energy 
per some time period is capped. One of the most 
common ways to define stability is bounded-input 
bounded-output (BIBO) stability. 


A system is said to be BIBO stable if, for any 
bounded input (meaning that the magnitude of the 
signal never exceeds some finite value), the output 
is also bounded (but not necessarily by the same 
value as the input). Mathematically, we can put it 
like this: 


Consider a discrete-time system H and arbitrary 
input signal x[n] with |x[n]|<M1€R, Vn. Let 

y[n] =H{x[n]}. H is BIBO stable if there exists some 
M2eR such that |y[n]| <M2, V n. 


BIBO stability is a guarantee, a stamp of approval on 
a discrete-time system, certifying that the output 
will always be capped by some value...so long as the 
input also is. 


The Impulse Response and BIBO Stability 
in LTI Systems 


For LTI systems, there is again a special relationship 
between the system and its impulse response. The 
system's impulse response can tell us two things 
regarding the system: if it is BIBO stable, and if so, 
the particular value that bounds the output. 


Consider a discrete-time LTI system H with impulse 
response h[n], and arbitrary input signal x[n] with | 
x[n]|<M1, V M. H is BIBO stable if and only if 
there exists some M2 such that |[h[n]|]1 = M2. 
Furthermore, this M2 (if it exists) bounds the output 
of the system: |H{x[n]}| <M1M2. 


Proof 


The proof of our "if and only if" statement has two 
parts. First, we must show that the existence of 
Ih[n]||1 implies BIBO stabilty. Second, we must 
show that BIBO stability implies the existence of 
Ih[n]||1. We'll start with the first part, the "if." Let H 
be a system with bounded impulse response 

Ih[n]||1 = M2. Consider the output corresponding to 
an arbitrary input x[n] which is bounded by M1: 


lyn] | = |H{x[n]}| = |2m = — c eh[n—m]x[m]| =| 
Ym = —  cx[n—mJh[m]| <2ym= — ~ o |x[n—m]]| 
h[m]|]<2m= — © oM1|h[m]]| <M1zm= — ~ | 
h[m]|<M1M2. 


The output is bounded by a finite value, M1M2. So, 
if the impulse response ||h[n]||1 for an LTI system 
exists, then the system is BIBO stable. 


The "only if" side of the proof is to show that if a 
system is BIBO stable, the norm ||h[n]|l1 of its 
impulse response must exist. We will demonstrate 
this by proving the contrapositive: if ||h[n]||1 is 
unbounded, the system is not BIBO stable. It will 
require the use of a the function sgn{x[n]}, which 
outputs 1 when x[n] is positive, — 1 when it is 
negative, and 0 when it is 0. 


Consider an arbitrary impulse response h[n] that is 
not absolutely summable, i.e., ||h[n]||1 is unbounded. 
For the system to be BIBO stable, any bounded input 
must result in a bounded output. So let the input 


x[n] =sgn{h[ —n]}. Clearly, x[n] is bounded: 
lIx[n]lleo =1. But 


ylO] =xm= — © oh[0O—m]x[m] =xym= 
— co ooh[O— m]sgn{h[m]} = Xm = — e  |h[m]| 
=|lh[n]|I1. 


So in this case y[0] is unbounded, so the system is 
not BIBO stable. 
Example of an impulse response of a BIBO system. 
Example of an impulse response of a system that is 
not BIBO stable. 


Examples of BIBO and non-BIBO Systems 


Suppose an LTI system (remember, impulse response 
absolute summability being equivalent to BIBO 
stability only applies to LTI systems!) has the impulse 
response h[n] = 1n2u[n—1]: 

h{n] 


As the norm ||h[n]||1 of this impulse response exists 
(Zn = 1|1n2|=7226), the system is BIBO stable. But 
suppose h[n] = 1nu[n—1]: 


For this system, the norm |lh[n]||1 does not exist, as 
Xn=1-|1n|—o. So this system is not BIBO stable. 


Because of this impulse response property, it is 
evident that all FIR systems are BIBO stable (for a 
finite sum of finite values is finite). 


Conclusion 


Once again we have seen the significance of the 
impulse response for LTI systems. They can be used 
to compute a system's output, and they can be used 
to determine the causality of a system. Here we 
have seen they can also be used to determine the 
BIBO stability of a system, a characteristic of great 
importance. But note two things. First, the 
connection between the impulse response and 
stability (as with those other connections) applies 
only to LTI systems. For example, the system 
H{x[n]} =1x[n]+1-—1 is not stable (the output 
explodes as the input tends to —1), yet its impulse 
response is —128[n], which is absolutely summable. 
And second, just because a system is not BIBO 
stable, this does not necessarily mean that it is not 


useful, but rather that care must be given to what is 
input into the system. 


Orthogonal Bases 


Recall from before, how we can understand discrete- 
time signals to be vectors in a vector space. There 
are some very useful reasons why we might want to 
express some signal in a vector space in terms of 
other signals in that space. To better understand 
how all of this works, and to give us some 
mathematical foundations for it, we will consider 
the concept of bases. 


The Basis of a Vector Space 


Suppose we have some vector space V, such as RN 
or CN, i.e., the set of all real or complex valued 
finite-length (of length N) discrete time signals. We 
define a basis for V as a set of vectors {bk},bkEV 
which span V and are linearly independent. By 
spanning V, we mean that any vector in V can be 
expressed as a linear combination of one or more 
vectors in {bk}k=ON—1: 


VxEV,x =2k=ON-—1lakbk = a0b0+a1b1+::-+aN 
—1bN-—1,akE€C. 


By the vectors in the set {bk} being linearly 
independent, we mean that no vector in that set can 
be expressed as a linear combination of any of the 
others. The number of these spanning and linearly 


independent vectors in the basis is the dimension of 
the basis. Our example bases RN or CN are of 
dimension--you guessed it--N. 


The Basis Matrix 


For the sake of cleaner and simpler mathematical 
expression--as well as the ability to connect the 
concept of a basis with other linear algebra tools-- 
we can create a matrix with basis vectors as 
columns. If the dimension of the basis is N, then this 
collection of basis vectors will be an N x N matrix, 
which we'll call B: 


B = [bO|b1|---|bN—1]. 


Recall we can express any vector in a vector space 
as a linear combination of the basis vectors. We can 
put these weights {ak} into an N x 1 column vector: 


a=[a0al:aN—1]. 


With the basis matrix B and the weights vector a, 
we can express the linear combination using a 
simple matrix multiplication: x =Ba. So we see that 
we can use the basis matrix and weights vector to 
refer to any vector in the vector space, through the 
linear combination of the basis vectors (which we 
can express with a matrix multiplication). Now, it is 
a natural question to ask, given some vector x in the 


vector space, how do we find the weights a that will 
produce the expression x = Ba? If we would like to 
express x as a linear combination of basis vectors, it 
is important that we know how to find those 
weights! Thankfully, it is very straightforward. 
Simply multiply each side of the matrix 
multiplication equation by the inverse of the basis 
matrix: 


x= BaB —- 1x = B— 1BaB— 1x=IaB-— 1x=a. 


The weight vector is simply the inverse of the basis 
matrix, times the vector x (the one we want to 
express in terms of the basis matrix). 


Orthogonal and Orthonormal: Special 
Bases 


We have seen that a basis is a special collection of 
vectors from some vector space: it is a collection 
that spans the space, and is mutually linearly 
independent. If we add another requirement or two, 
we end up with two important sub-classes of bases. 
Suppose the vectors in some basis are not only 
spanning and linearly independent (which of course 
they must be, by definition, to form a basis), but 
that they are also mutually orthogonal, that is, that 
the inner product of any basis vector with any other 
basis vector is 0: 


V bke{bk}k=ON —1,(bk,bl]) = 0,k 1. 


If such is the case, then this basis is said to be an 
orthogonal basis. 


So an orthogonal basis is a particular kind of basis, 
one whose vectors are mutually orthogonal. Among 
orthogonal bases, there are some whose vectors all 
have unit 2-norms: 


(bk,bl) = 0,k #I|[bkII2 = 1. 


A basis with this additional property is known as 
an orthonormal basis. 


Basis Matrix of an Orthonormal Basis 


Like any other basis, the vectors of an orthonormal 
basis can be put together to form a basis matrix. 
Recall how we find the weights a to express some 
vector x in terms of the basis matrix B: a=B— 1x. 
The reason that orthonormal bases are so special is 
that, in contrast to other bases, their matrix inverses 
are extremely easy to find. For orthonormal basis 
matrices, B— 1=BH. The inverse of an orthonormal 
basis matrix is simply its Hermitian (conjugate) 
transpose! So for these bases, a= BHx. In linear 
algebra, a matrix that has the property that its 
inverse is simply its Hermitian transpose is called a 
unitary matrix. If the matrix is real-valued, it is 


called an orthogonal matrix (this is a bit of 
unfortunate nomenclature, considering its columns 
are actually mutually orthonormal), and its 
Hermitian transpose is simply the regular matrix 
transpose, BT (sometimes written as B’). 


Orthonormal Basis Signal Representation 


Putting this all together, we can see the two key 
aspects of signal representation with orthonormal 
bases. There is the synthesis side, that we can build 
up any vector x in the vector space through a linear 
combination of basis vectors. And there is the 
analysis side, that we can find the proper weights 
of this linear combination by simply multiplying x 
by the Hermitian transpose of the basis matrix. 


* Synthesis: x = Ba = Xk=ON-— 1akbk 
* Analysis: a= BHx ,or, ak = (x,bk) 


Eigenanalysis of LTI Systems (Finite-Length Signals) 


In the study of discrete-time signals and systems, 
concepts from linear algebra often provide 
additional insight. When it comes to LTI systems, a 
certain area of linear algebra is particularly helpful: 
eigenanalysis. 


Eigenvectors and Eigenvalues 


Given a square matrix (one that has the same 
number of rows as columns) A, a vector v is an 
eigenvector with corresponding scalar eigenvalue 7 
if: Av=Av. There is a geometric interpretation to 
this eigenanalysis of the matrix A. Multiplying a 
matrix by one of its eigenvectors produces simply a 
scaled version of that same eigenvector (scaled by a 
factor of X.), so a matrix multiplication of an 
eigenvector does not change its orientation, only its 
strength. 


Consider this example in two dimensions, a square 
matrix A: 


A=[3113], 
and the vector v: 


v=[1-1]. 


Note what happens when we multiply the matrix A 
by the vector v: 


Av=[3113][1—-1]=[(3—1) 
(1 —3)] =[2-—2] =2[1—-1] =2v. 


So running the vector v through the matrix A simply 
scales the vector by 2. 


Eigendecomposition: Handling Multiple 
Eigenvectors and Eigenvalues 


We have seen what eigvenvectors and eigenvalues 
are for square matrices. Now, an N XN matrix will 
have N eigenvectors (not necessarily distinct), each 
with its own eigenvalue. We can put all of these 
vectors and values into their own matrices. Suppose 
the eigenvectors of the matrix are {vm}m=ON-—1 
and the values are {Am}m=ON-— 1. Then we can 
organize all of them like this: 


V=[v0|v1|---|vN—1]A=[AOA1.AN-1]. 


With those vectors and values collected like that, we 
can express the eigenvector/value property Av=AVv 
for all of the eigenvectors and eigvenvalues of the 
matrix A at once: AV=VA. 


Diagonalization 


Now, if the square matrix V is invertible (which will 
be the case if the columns of A are linearly 
independent, which of course happens if the vectors 
form a basis), then we can do some special things 
with it. 


Recall the eigendecomposition 

relationship: AV=VA. If V is invertible, then we can 
multiply each side of the equation by V—1: V 
—1AV=V—-1VA=IA=A. So V—1AV=A. Because 
multiplying A by V on one side and V—1 on the 
other produces a diagonal matrix, we say that the 
matrix V diagonalizes A. If we multiply the 
eigendecomposition matrix the other way, we will 
have that A=VAV—1. One of the reasons (more of 
which we'll see later) why we consider this 
diagonalization of A is that it is easier to matrix 
multiply with diagonal matrices than full ones. 
When a complex harmonic sinusoid is input into an 
LTI system, the output is a scaled version of the 
input. Here the real (cosine) and imaginary (sine) 
parts of the sinusoid are plotted as the input. Note 
how the system merely scales the inputs. 


LTI Systems and Eigenanalysis 


Perhaps you may be wondering how all of this 
linear algebra relates to discrete-time systems. For 
LTI systems operating on finite-length discrete-time 


systems, the input output relationship is: y = Hx, 
where H is a circulant matrix (each row being a 
circularly shifted version of the system impulse 
response h). 


Let's see what the eigenvectors of H are. These will 
be the finite-length signals that, when input into the 
system, emerge as outputs simply as scaled versions 
of themselves. In that sense, they are somehow 
fundamentally related to all LTI systems. 


It so happens that, remarkably, any and all LTI 
systems for finite-length (length N) signals have the 
exact same set of eigenvectors! The eigenvectors for 
any LTI length-N system are complex harmonic 
sinusoids: 


sk[n] = ej2aNknN = 1N(cos(2xNkn) + jsin(22Nkn)),0 <n. 
—1. 


So, if we have an LTI system--any LTI system--then 
giving an sk as an input will result in the output 
being AKsk, with the particular values of Ak of 
course being dependent on the system: 


A, cos 2Zkn), k = 2 


a aS ae At; | 
to i]t 10 15 aI 
| Ag sin(22kn), i 


To prove this special property of LTI systems, we 
simply compute the circular convolution sum for an 
LTI system with arbitrary impulse response h[n] and 
input of the form ej2mNknN: 


sk[n]@®h[n] = Xm =ON — 1sk[(n 

—m)N]h[m] = Xm = ON — lej2aNk(n 

—m)NNh[m] = xm =ON — lej2xNk(n 

—m)Nh[m] = »m=ON — lej2xNknNe 

—j2nNkmh[m] =(%m=0ON-— le 
—j2nNkmh[m])ej2xNknN= Aksk[n] , AK=(2m=0N 
—1le—j2nNkmh[m]). 


This proof reveals how we are to find the 
eigenvalues (Ak) that correspond to each harmonic 
sinusoid eigenvector: they are simply the inner 
products of the eigenvectors with the system's 
impulse response. Each value Ak is called the 
system's frequency response at frequency k, because it 
indicates how the system scales inputs of that 
particular frequency. It is a significant enough 
characteristic of the system to warrant its own 


notation: H[k]. 

Graphical representation of the real and imaginary 
parts of the discrete-time finite length LTI system 
eigenvector matrix S.The real part of the 
eigenvector matrix S: cos(2sNkn)/N. The real part 
of the eigenvector matrix S: sin(2xNkn)/N. 
Eigendecomposition of discrete-time finite length 
LTI systems. 


Eigendecomposition of LTI Systems 


As with matrices in general, we can apply an 
eigendecomposition on an LTI system matrix. For 
LTI finite-length systems, these matrices are 
circulant: 


y H Ny 
LT I system 


We have seen that the eigenvectors of LTI systems 


are harmonic complex sinusoids. We can stack these 
up into a single matrix S, the entries of which are 
Sn,k, which is plotted below: 


0 


k 


Likewise, we can plot the respective eigenvalues of 
the eigenvectors, which above we defined to be the 
values of the frequency response of the system: 


A=[AOA1:.AN—1] = [Hu[0]Hu[1]°--Hu[N—-1]]. 


Putting the equation y = Hx together with the 
sor gaa H= Bast, we have: 


We already know one way of understanding how 
LTI systems operate on signal inputs: they convolve 
them with the system's impulse response 
(represented in linear algebra form by the equation 
y = Hx, where H is circulant). The 
eigendecomposition gives us another understanding. 
The matrix SH takes the input and extracts what 
would be the coefficients of the input's 
representation as a linear combination of harmonic 
sinusoids (it turns out this is called the discrete 
Fourier transform). Then, multiplication by the 
diagonal matrix A modifies these coefficients in a 
way that is particular to the system H (all LTI 
systems have the same S and SH). Finally, 
multiplication with the matrix S takes the modified 
coefficients and expresses them as a linear 


combination of harmonic sinusoids to give the 
output y (it turns out that operation is called the 
inverse discrete Fourier transform). 


The Discrete Fourier Transform 


The Significance of Complex Harmonic 
Sinusoids 


Recall that complex harmonic sinusoids have the 
form: 


sk[n] = ej2aNknN. 


In the context of finite-length discrete-time signals 
and systems, complex harmonic sinusoids are 
special, for several reasons. First, the collection of 
them {sk}k=ON-—1 form a basis for CN, meaning 
that any signal in CN can be represented as a linear 
combination of that set of N vectors. Not only this, 
but they are an orthonormal basis, which means 
that the weights of the linear combination can be 
found simply through the inner product of the signal 
with the harmonic sinusoid in question. And finally, 
these signals {sk}k =ON—1 are important because 
they are eigenvectors of finite-length discrete-time 
LTI systems. This means that whenever a complex 
harmonic sinusoid is given as an input to any LTI 
system, the output is simply a scaled version of that 
input. 


All of these special features come together as we 
now consider a very important component of the 


study of signals and systems. 
A discrete-time finite length signal x[n], and the 
magnitude of its DFT, |X[k]]. 


The Discrete Fourier Transform 


In the early 19th century, Jean Baptiste Joseph 
Fourier showed that any function (later 
mathematicians would more rigorously qualify that 
statement) could be composed as a (possibly 
infinite) sum of harmonic sinusoids. This work 
resulted in what would be an entire branch of 
mathematics, Fourier analysis. Fourier analysis 
extends to many different kinds of signals, including 
discrete-time finite-length signals. For those, the 
analysis produces the discrete Fourier transform 
(DFT). For signals x[n]€CN, the normalized DFT 
and inverse DFT are: 


X[k] = xn = ON — 1x[nJle —j2aNknNx[n] = xk=0ON 
—1X[k]ej2xNknN. 


We see that taking the inner product of x[n] with 
harmonic sinusoids of different frequencies k (which 
is what the first sum represents) produces a series of 
frequency coefficients X[k]. That is the DFT. We can 
also say that it is the analysis aspect of the Fourier 
transform, for it gives us a frequency analysis/ 
breakdown of the signal x[n]. The frequency 
coefficients can be used to reconstruct x[n] using 
the second sum, which is the inverse DFT. It is 


known as synthesis, for it shows how x[n] can be 
built up as a linear combination of harmonic 
complex sinusoids, with the coefficients X[k] telling 
us how much of each one is needed. 


x{n] 
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We can express the above signal analysis and 
synthesis in matrix notation. First we stack up the 
complex harmonic sinusoids {sk}k =ON—1 as 
columns in a matrix S: 


S=[s0|s1|---|sN—1]. 


With these vectors in a matrix, then the signal x is 
composed of the linear combination of the sinusoids 
(the synthesis operation), with the weights in a 
vector X, through the matrix multiplication: x = SX 
(synthesis; inverse DFT). Finding the weights X, 
given x, requires the inverse of the basis matrix S. 
However, since S is unitary, the inverse is simply 
the Hermitian transpose: X = SHx (analysis; DFT) 


Normalized and Un-normalized DFT 


To this point, we have been working with the 
normalized expression of the DFT and its inverse: 


¢ DFT: X[k] =xn=ON — 1x[n]e—j2aNknN 
¢ Inverse DFT: x[n] = xk =ON — 1X[k]ej2aNknN 


This form of the DFT has a certain symmetric 
elegance to it, the only difference in the DFT and its 
inverse being the negation of the exponent in the 
complex harmonic sinusoid. 


However, it is far more common to use an un- 
normalized definition of the DFT, one which puts a 
IN scaling factor on the inverse, rather than the 
normalized version splitting this with a 1N in each 
formula. In practice, the DFT is virtually always 
understood in terms of the following un-normalized 
definition: 


¢ DFT: X[k] =xn=ON —1x[n]Je—j2aNkn 
* Inverse DFT: x[n] = 1NXk=ON — 1X[k]ej2aNkn 


Interpretations 


So what exactly does the DFT mean? There are a 
variety of interpretations for what the DFT does. We 
have already explained it in terms of analysis and 


synthesis. With the inverse DFT, we can build up 
any signal x[n] as a weighted sum (also known as a 
linear combination) of complex harmonic sinusoids; 
this is known as synthesis, which is in accordance 
with the Merriam-Webster dictionary definition of 
synthesis as being "the composition or combination 
of parts or elements so as to form a whole." But how 
are we to know how much of each sinusoid goes in 
to the synthesis? The DFT tells us the coefficient 
values, one for each of the N sinousids. This is the 
analysis aspect of the DFT, for it "analyzes" the 
signal in terms of the frequencies in it, how much 
there is of every possible frequency within the 
signal. To use a cooking metaphor, the analysis tells 
us the amount of each ingredient that goes into the 
dish that is our signal, and the cooking process is 
the synthesis which takes the ingredients and 
creates the finished product. 


There is another way of looking at the DFT and 
inverse DFT. We note that there is a one-to-one 
correspondence between any signal x[n] and its DFT 
X[k]; a signal x[n] has only one DFT, and any X[k] 
has only one inverse. So x[n] and X[k] are really 
referring to the same one thing--some signal--but 
only in different ways. x[n] describes the signal in 
terms of its value at every given time location n, 
while X[k] describes the signal in terms of how 
much of each frequency contributes to it. So there is 
a time-domain description of the signal, which is 
x[n], and a frequency-domain description of it as 


well, which is X[k]. We could think of the time- 
domain and frequency-domain as two different 
languages, like English and Spanish. For some kind 
of idea or expression, it can be represented in either 
one language or another, and given the expression 
in one language, it can be translated to the other. Of 
course, this is not a perfect metaphor; spoken 
languages are not always one-to-one (and some 
words in one language do not have an adequate 
word in another language). What this interpretation 
does convey is that time and frequency are two 
different, and yet also equivalent, ways of 
expressing a signal entity. 

Heideman, Michael, Don Johnson, and C. Burrus. 
"Gauss and the history of the fast Fourier transform." 
IEEE ASSP Magazine 1, no. 4 (1984): 14-21. The a) 
crank, b) gears, c) arms, d) levers, and e) springs of 
a harmonic analyzer can together either synthesize a 
function or analyze it, writing the resultant plot on a 
f) piece of paper. 


Excursus: A Physical Implementation of 
Fourier Analysis and Synthesis 


The inverse DFT shows us that a time-domain signal 
can be represented as a summation of complex 
sinusoidal functions, while the DFT reveals how 
much of each complex sinusoid is required in the 
summation. Before the advent of computers, 
calculating either the DFT or inverse DFT was 
obviously a tedious process. Even with the use of 


some algorithmic shortcuts (e.g., by Carl Friedrich 
Gauss in the early nineteenth century[footnote]), 
the calculations had to be done by hand. 


But "by hand" doesn't have to mean writing out the 
equations and consulting mathematical tables. In 
the early twentieth century, Albert Michelson 
invented a machine that could perform Fourier 
synthesis and analysis. This "harmonic analyzer" 
used a series of gears, levers, and springs to 
physically implement the inner products we see in 
the DFT and inverse DFT equations. 
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Recall that in synthesis, the inverse DFT adds 
together a weighted sum of complex sinusoidal 
functions. The machine performs a similar 
operation, with each part of the machine 
corresponding to an aspect of the inverse DFT 
formula. The crank of the machine rotates a series of 
gears that each turn another gear. Each of these 


gears would correspond with a certain k of the 
formula, with the turning of the crank being the 
equivalent of progressing through the sinusoids in 
time (the variable n). The rightmost gear would 
have a single revolution per turn of the crank, then 
the next one over, two revolutions, and so on. All of 
these gears are connected with rods to a set of arms. 
As the gears turn, the arms rock up and down. Thus 
turning the crank produces a set of sinusoids of 
varying frequency (in the figure above, from k=1 to 
k=80!). 


Of course, the synthesis operation must weight these 
sinusoids. The sinusoids are weighted (the X[k] 
values in the inverse DFT formula) by setting each 
lever at a specified position along the rocker arm. 
The farther towards either the edge of the arm, the 
more the lever will move along with the arm. At the 
top of the machine, the ends of these levers are 
attached to springs that are then attached to a bar 
with a single counter-spring. All of these springs 
being connected to the same pivoting bar "sums" the 
movement of all of the levers, which corresponds to 
the sum ~ in the formula. At the end of the pivoting 
bar is a lever that is tied to a pen, which moves up 
and down as a piece of paper tranverses left to right 
while the crank turns. Thus the pen produces on the 
paper the output of the inverse DFT summation. 


The machine can also be used to perform analysis. 
In order to find the proper weights that would 


produce a desired function through analysis, a single 
period of the function would be "inputted" into the 
machine by setting the lever positions across all of 
the arms into the shape of the function. That shape 
across all the arms corresponds to the x[n] of the 
DFT formula, and now the number of arms 
correspond to the time variable n. As the crank is 
turned, the arms will make a shape of a sinusoid (as 
you can see in section c) of the figure). One turn of 
the crank will produce a very long wave, a single 
sinusoidal period across all the bars, i.e. a sinusoid 
with frequency k=1. The shape of the function 
weighted across all of these bars corresponds to the 
weighting in the DFT formula, and once again the 
pivoting bar will sum everything together. After one 
turn, the position of the pen will correspond to 
X[1], after two turns, X[2], and so on, thereby 
computing the DFT. 


To see a video of the harmonic analyzer at work, 
visit http://www.engineerguy.com/fourier/. 


Discrete Fourier Transform Properties 


Recall, for discrete-time finite length signals, the 
definition of the DFT and the inverse DFT (given 
here in its normalized form): 


¢ DFT: X[k] =xn=ON — 1x[n]e—j2aNknN 
¢ Inverse DFT: x[n] = Xk =ON — 1X[k]ej2aNknN 


Recall also its much more commonly used un- 
normalized form: 


¢ DFT: X[k] =xn=ON —1x[n]Je—j2aNkn 
* Inverse DFT: x[n] = 1NXk=ON — 1X[k]ej2aNkn 


A signal x[n] and its DFT X[k] (recall there is a one- 
to-one correspondence) are referred to as a DFT 
pair. The DFT has a variety of properties, which we 
will now consider. 

The DFT X[k] of any N-length signal is periodic, 
with period N. 


The DFT and Its Inverse Are Periodic 


The DFT is defined for finite-length (length N) 
signals; so, for x[n], n runs from 0 to N—1, as does 
k. But let's see what happens when we consider a 
value of k outside this range in the DFT formula, say 
X[k+IN], where | is some nonzero integer: 


X[k+1IN]= X%n=ON — 1x[nJe —j2xN(k 
+IN)n=%n=0ON — 1x[n]Je—j2aNkne 

—j2nNINn = Xn =ON — 1x[nJe —j2aNkn-1= Xn=O0N 
— 1x[n]e —j2xNkn=X[k]. 


As X[k + IN] =X[k], the DFT is periodic: 
|X [| 


By the same token, the inverse DFT is also periodic: 


x[n+IN] = INXk=ON — 1X[k]Jej2xNk(n 
+ IN) =%k = ON — 1X[k]Jej2aNknej2aNkIN = Xk = ON 
— 1X[k]ej2aNkn-1 = Xk =ON —1X[k]ej2aNkn=x[n]. 


Again this is to be expected because the complex 
harmonic sinusoids of the inverse DFT sum are 
periodic. This also further illustrates the fact that 
any finite-length signal can be understood as a 
single period of a corresponding periodic signal. 

For an N= 16 length signal, one way to express the 
range of frequencies is for k to run from 0 to N—1, 
which corresponds to frequencies between 0 and 2z. 
For an N= 16 length signal, another way to express 
the range of frequencies is for k to run from —N2 to 
N2-—1, which corresponds to frequencies between 
—s and x. 


DFT Frequency Ranges 


A complex harmonic sinuoid ej(2sNk)n has, by 
definition, a frequency of 2xNk, which may label as 
wk. In the DFT of a signal of length N, the variable k 
ranges from 0 to N—1, values that correspond to 


frequencies from 0 to (just about) 27: 
|X [k}| 


0 5 10 15 


However, as we saw above, since the DFT is 
periodic, 0 to N—1 is not the only range of k we 
may use to express the DFT. Since X[k] =X[k—N], 
we may let k run from —N2 to N2—1 (for even N, 


that is; for odd N it would be -N—12 to N—12): 
|X[k]| 


Shifts in Time and Frequency 


Let x[n] and X[k] be a DFT pair (i.e., X[k] is the 
DFT of x[n]). A circular shift in time on x[n] will 
produce a phase shift on X[k]: 


x[(n —m)N] «?DFT e—j2aNkmX[k]. 


To prove this relationship, we first note that for the 
circular shift (n—m)N, there is some integer | such 
that(n—-m)N=n—m+IN. We will use that fact for a 
change of variables in our proof: 


DFT{x[(n —m)N]}=*n=ON-1x[(n—m)N]e 
—j2nNknLet r=(n—m)N=n—m+IN=2r=0N 

— 1x[rJe —j2aNk(r +m—IN)=2r=ON-—1xI[rle 
—j2aNkre —j2mNme — j2NIN =e —j2aNm>2r = 0N 
— 1x[r]Je —j2nNkr-1 =e —j2nNkmX[k]. 


As you might expect from the symmetrical 
similarities between the DFT and inverse DFT, a 
comparable relationship exists in the other direction 
as well. A circular shift in the frequency domain of a 
signal results in the modulation of the signal in the 
time domain: 


ej2NInx[n] <?DFT X[(k—I)N]. The proof of this 
property follows the same approach as that of the 
first. 


The DFT is Linear 


As the DFT of a signal is essentially a weighted sum, 
it follows naturally that the DFT operation is linear. 

So if we have two arbitrary signals x1[n] and x2[n] 

(with DFTs X1[k] and X2[k]) and arbitrary 


constants al and a2, then the DFT of 
a1x1[n] +a2x2[n] is simply a1X1[k] + a2X2[k]: 


x1[n] <°DFT X1[k], x2[n] «?DFT X2[k]alx1[n] +a2x[2] 
The proof is straightforward: 


DFT{alxl[n]+ta2x2[n]} =2n=O0ON 
—l(alxl[n] +a2x[2])e-j2nNkn 
=(altn=ON-1xl1l[nle-j2xNkn) 
+Ca2in = ON —)x2[nhe-=—j2aNkn) 
= 01X11. hee x2] 


The Convolution/Multiplication Time/ 
Frequency Relationship 


We now reach perhaps the most significant DFT 
property, the relationship between convolution and 
multiplication in time and frequency. Suppose we 
have two discrete-time finite length signals x1[n] 
and x2[n], whose DFTs are X1[k] and X2[k]. Let 
y[n] be the circular convolution 

y[n] =x1[n]@x2[n]. Then the DFT of y[n] is 
equivalent to the product of the DFTs of x1[n] and 
x2[n]: Y[k] =X1[k]X2[k]. Or, to express this 
relationship in DFT pair notation, we have: 


xl{(n]@®x2[n]-DFTX1[k]X2[k] 


The proof of this relationship will use the r=(n 


—m)N=n-—m+IN change of variable we used 
earlier: 


DFT{x1[n]@x2[n]} = Xn=ON-—1(2m=ON-1x1[(n 
—m)N]x2[m])e —j2aNkn = Xm =ON — 1x2[m] 

(Xn = ON — 1x1[(n—m)N]e—j2nNkn) = Xm =0N 
—1x2[m](®r=O0N — 1x1[r]e —j2aNk(r+m 

+IN)) =xm =ON — 1x2[m](“r= ON — 1x1 [rle 
—j2aNkre —j2nNkme — j2NkIN) = (2m = ON 
—1x2[m]e—j2nNkm)(“r=0ON—1x1[r]e 

—j2maNkr) = X1[k]X2[k]. 


So we see that convolution in two signals’ time 
corresponds to their multiplication in frequency. It 
is certainly an interesting relationship, but it is also 
very practical: N multiplications in the frequency 
domain are obviously fewer than the N2 
multiplications required to compute for an N-length 
signal. Of course, this requires first that the DFT of 
the two signals be computed (and then the inverse 
DFT of the product), but we will see there are ways 
to compute that efficiently. What it all means is that 
a system's output (which is found via convolution) 
can be computed efficiently by way of 
multiplication in the frequency domain. 

A real-valued finite-length signal x[n], with the real 
and imaginary parts of its DFT.For real-valued 
signals such as x[n] here, the real part of their DFTs 
are even and the odd part is odd, which is more 
readily seen when they are plotted with k ranging 
from —N2 to N2. 


DFT Symmetry Properties 


As the DFT is a weighted sum of complex signals 

(complex harmonic sinusoids), it follows that the 

DFT of signals is, in general, complex-valued. 
x{n] 
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Perhaps you have noticed, in the figures above, that 
there is a certain symmetry for these (real-valued) 
signals' DFTs. For the real part of the DFT, you see 
that Re [X[k]] =Re [X[N-—k]], and that Im [X[k]] = 
—Im [X[N—k]]. Recalling that the DFT is periodic, 
if we consider frequencies from — N2 to N2, then we 
have Re [X[k]] =Re [X[—k]], and that Im [X[k]] = 
—Im [X[—k]]. Or, in words: the real part of the 
DFT is even, and the imaginary part is odd. 


x(n] 
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This symmetry was not a coincidence for the DFT of 
that particular x[n]; because the complex harmonic 
sinusoid signals that make up the DFT sum formula 
are conjugate symmetric (the real part is even, the 
imaginary part is odd), the DFTs of all purely real 
signals will also be conjugate symmetric. The 
converse is true for purely imaginary signals; their 
DFTs will be symmetric such that the real part is 
odd and the imaginary part is even. If, in addition to 
being purely real or purely imaginary, the signal 
itself is also even or odd, that will further limit the 
characteristics of its DFT. Note that if a signal is 
complex-valued, its DFT will not--in general--exhibit 
any such symmetries. 
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The Fast Fourier Transform 


Tasks and Algorithms 


The first thing to know about the fast Fourier 
transform, or FFT, is that it actually is not a unique 
type of mathematical transform at all. The FFT is 
simply a method, or algorithm, for the task of 
computing a DFT. In that way, it is similar to the 
task of sorting a group of numbers from smallest to 
largest. No matter how you might go about the task 
of sorting the numbers, there is only one correctly 
sorted answer. Some methods to sorting, however, 
may be quicker (take fewer steps) than others. For 
example, you could put the numbers in one list, and 
then progressively traverse through the list, 
comparing two adjacent items at a time and 
swapping them if necessary; each time you traverse 
the list, you will have put the next largest number in 
its place. This approach is usually called a "Bubble 
Sort." Another approach is to split the list of 
numbers into two sets, then divide each set in two, 
and divide those sets, and so on, until each little set 
has two (or fewer) numbers. You will then sort each 
set of two numbers (easy!), and then progressively 
combine these sorted sets together (also easy!) until 
a single sorted set remains. This is called a "Merge 
Sort," and is known as a "divide and conquer" 
algorithm. 


Now, these two sorting methods achieve the same 
results, but the Merge Sort typically does so in far 
fewer steps (although the Bubble Sort will process 
an already or nearly sorted list quicker). If there are 
N items in a list, Merge Sort requires about Nlog2N 
steps, while the Bubble Sort requires about N2. Even 
with lists as small as 32 items, that makes a big 
difference, the difference between 160 and 1024. So 
it is with the discrete Fourier transform. The DFT is 
a simple formula: 


X[k] = =n = ON — 1x[n]Je —j2aNkn. 


It would be possible to calculate this as a sum of N 
multiplications (going from n=0 to n=N-—1), doing 
that N times (for each k, again from 0 to N—1). 
That's about N2 multiplications and additions. But, 
there are many different ways one could calculate 
the DFT, just as there are many ways to sort a list. 
As with Merge Sort, there is a divide-and-conquer 
approach to finding the DFT which does so with far 
fewer computational steps, only about Nlog2N. 

The FFT is a "divide and conquer" algorithmic 
approach to finding a DFT. An N-length DFT can be 
found as the weighted sum of two N/2-length DFTs, 
and each of those are the sum of two N/4-length 
DFTs, and so on.The DFT of a length-8 signal can be 
split into the the weighted sum of two length-4 
DFTs, namely, the DFTs of the even and odd indices. 
Each length-4 DFT is split into the weighted sum of 
two length-2 DFTs. The length-2 DFTs are found, 


and from there the rest of the overall DFT can be 
computed according to the half-size weighted sums. 
The 2-length DFTs are simple to find, just an 
addition and a subtraction. When those operations 
are represented graphically, they are said to have a 
"butterfly" structure. 


The FFT: A Divide and Conquer Algorithm 
for Finding a DFT 


The intuition behind the FFT algorithm is the same 
as that of the Merge Sort. The Merge Sort operates 
on the idea that it is easy to take two sorted lists 
and combine them into one large sorted list: you 
simply compare the beginning of each list, pick off 
the smallest of the two and place it in a new list, 
and then keep on doing that until the two lists are 
empty. 


Something similar happens with the DFT. Suppose 
you have some signal x[n], which you can split into 
two smaller signals in terms of its even and odd 
indices, x[2n] and x[2n+ 1]. Suppose now that you 
already know the DFTs of these two smaller signals, 
and call them E[k] and O[k]. Then it is easy to find 
the DFT of the whole signal, X[k], from these 
smaller ones: 


X[k] =E[k] + ej2aNkO[k]. 


A proof of this follows. For simplification of 


notation, we let WN =e-—j2aN: 


X[k] =xn=0ON 

— 1x[n]WNkn = Xn =0N/2 — 1x[2n]WNk(2n) + Xn=ON/2 
+1]WNk(2n 

+1)=XYn=O0N/2-1x[2n]WN2kn + WNk2Yn=O0N/2—-1x[2 
+1]WN2kn = Yn =0N/2 —1x[2n]WN/2kn + WNkxn=ON, 
+1]WN/2kn=E[k] + WNKO[k]. 


You will note that the k in X[k] runs from 0 to N 
—1, but that the even and odd DFTs are only of 
length N/2. It may seem a little strange to consider 
the values of these DFTs for k greater than N/2 (as 
is necessary in the formula). But there is no problem 
with this, since DFTs are periodic; as those DFTs are 
length N/2, we have E[N/2+k] =E[k] for those 
values between N/2 and N —1. So, we have that the 
DFT of an N-length signal is simply the weighted 
sum of two N/2-length DFTs. If we, somehow, are 
given the DFT of those even/odd sub-signals, there 
are only N additions and multiplications involved in 
finding the DFT of the whole signal, which is of 
course much better than the N2 needed to find the 
DFT by straightforward calculation of the DFT sum. 


Of course, this requires that we had the DFTs of the 
sub-divided signals! How do we find those? Well, 
again we will divide those sub-signals by even/odd 
indices, and then do that again, and again, and 
again, until we have lots of length-2 sub-signals. The 
DFT of a length-2 signal is very easy to compute. If 


we call such a signal t[n], then TLO] =t[0]+t[1] 
and T[1] =t[0] —t[1] (this is easy to verify for 
yourself). 


With the DFTs of our many length-2 signals in hand, 
we use the even/odd sum formula to combine them 
find the DFT of the length-4 signals, and then 
length-8, and so on until we have the DFT of the 
whole signal. This is known as a decimation-in-time 
FFT, and while it assumes a signal length that is a 
power of 2, there are ways to apply the approach to 
signals of other lengths. The figure below illustrates 
the process of splitting up the DFT into smaller and 
smaller half-size DFTs: 
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Computational Savings of the FFT 


Recall that a straightforward sum and multiplication 
computation of the DFT requires about N2 
operations; to be precise, it is N2 multiplications 
and N2—N additions. The FFT approach requires 
about Nlog2N operations: Nlog2N —N 
multiplications and Nlog2N additions. For small 
signal sizes, the difference between N2 and Nlog2N 
is not that significant. But as N gets even moderately 
large (say, N= 32 and up), the difference becomes 
larger and larger. By the time we consider typical 
signal sizes in signal processing, in which N is in the 
millions, taking N2 complex operations is 
prohibitive in terms of time and memory costs, 
while Nlog2N is very manageable. 


Fast Convolution with the FFT 


We have seen, previously, these four important 
signal processing facts: 


* the output of infinite-length LTI systems can be 


found via linear convolution 

¢ linear convolution can be found through 
circular convolution, through zero-padding 

¢ there is an equivalence between circular 
convolution in the signals' time domain and 
multiplication in their frequency domain 

¢ the FFT is an algorithm that can quickly 
compute a DFT 


Perhaps at the time it may have seemed that these 
findings were unrelated to each other. But as we 


now string them together, you can see the incredibly 


significant consequence: the output of LTI systems 
can be computed very efficiently. This truth has 
supported the incredible advances in signal 
processing over the past fifty years. 

Through zero-padding and forward and inverse 
DFTs, using the FFT algorithm, an LTI system's 
output can be found efficiently. 


Putting It All Together... 


Suppose we have an LTI system with an impulse 
response of h[n] and an input of x[n]. We would 


like to find the output, y[n]. We could find this 
output via linear (infinite length) convolution. If the 
length of the impulse response is Nh and that of the 
input signal is Nx, then about NhNx operations 
would be required to compute this convolution. It 
would also be possible to compute the output 
through circular convolution, by zero-padding each 
signal to be of length Nh+Nx—1, and then 
circularly convolving the zero-padded signals. This 
would not save any computational steps, but it is a 
significant insight because the circular convolution 
could be performed by using DFTs: simply take the 
DFT of each zero-padded signal, multiply the two 
DFTs together, then take the inverse DFT of the 
result. Again, this does not at first seem to save any 
computational steps, and in fact seems to add even 
more, except that we have seen that the FFT is able 
to perform DFTs in about Nlog2N operations. What 
that means is that if we zero-pad x[n] and h[n], 
then take the DFT of each (using the FFT 
algorithm), multiply these two together, then take 
the inverse DFT of the result, we can find the system 
output in about 2(Nh + Nx— 1)log2(Nh+ Nx-— 1) 
operations. As signal lengths increase, this can end 
up being huge computational savings over the NhNx 
operations it would take to find the output through 
a convolution sum. Below is a graphical depiction of 
the steps to quickly find the output, with the help of 
the FFT: 


yn 


Other Orthogonal Bases 


We have seen that the discrete Fourier transform 
(DFT) is an example of an orthogonal basis for finite 
length discrete-time signals. As a basis, a linear 
combination of its basis elements {bk} can represent 
any signal x: 


x = Yk =ON— lakbk =Ba. 


Furthermore, as an orthogonal basis, its basis 
elements are mutually orthogonal, meaning that the 
weights needed to build up a particular signal x can 
be found with a weighted inner product of x and the 
basis element in question. If all of the basis elements 
are of unit-norm, then the basis is orthonormal, and 
a simple inner product will find the weights: 


ak = (x,bk) , a= BHx. 


The DFT is a particularly useful basis, for its basis 
vectors ej2stNkn are the eigenfunctions of all LTI 
systems. As such, they are simply scaled when 
passing through a system. But as useful as the DFT 
is, there are other significant orthogonal bases to 
consider, as well. 

The real part of the DFT functions. The imaginary 
part of the DFT functions. The DCT functions are 
real-valued. For the this particular signal x[n], note 
now the real and imaginary parts of its DFT have 
significantly more large values than its DCT does. 


The Discrete Cosine Transform 


Since the DFT uses complex harmonic sinusoids to 
analyze and synthesize signals, the DFT coefficients 
of real signals are almost always complex-valued 
(unless the signal is even, in which case the 
coefficients are purely real). This normally poses no 
problem, but there are some circumstances in which 
real coefficients are preferable. One example is 
when signal compression is desired. In such 
applications, a signal is transformed and only a few 
of the largest transform coefficients are kept. Having 
to store both the real and the imaginary parts of the 
coefficients would mean twice the storage. It is for 
this reason that the discrete cosine transform (DCT) 
is used in signal compression applications (such as 
the JPEG image compression standard). The DCT is 
very similar to the DFT, as it also uses sinusoids as 
basis elements, except that its basis elements are 
purely real-valued. There are several varieties of the 
DCT, but the most popular is the DCT-II, which has 
the following basis elements: 


dk[n] = cos[aN(n+12)k],0<n,k<=N-1. 
d;(n| 
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Unlike the DFT basis functions, the DCT functions 


are completely real-valued. 
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Since the DCT basis functions are real valued, the 
DCT transform coefficients for a real signal are also 
real valued. As a result, there are cases in which 
signals can be represented reasonably well when 
reconstructed with just a handful of DCT 
coefficients, compared to the same number of DFT 
coefficients. 
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Note have the wavelet basis signals above have local 
sinusoidal features. They are thus useful for 


representing signals with similar characteristics. 


Wavelet Transforms 


We have seen that with the DFT (and DCT as well), 
signals can be represented in either the time 

domain (e.g., as x[n]), or in the frequency domain 
through their DFT (e.g., X[k]). Both representations 
are equivalent, of course, for one can always be 
derived from the other. But sometimes one way of 
expressing a class of signals is desired, such as when 
as signal can be expressed with just a few values. 
For example, the signal x[n] =5[n] can be expressed 
with just a single nonzero value in the time domain, 
but it does not have any nonzero values in the 
frequency domain. Thus we say the signal is 
localized in time. Likewise, the 

signal x[n] = ej2716n has only a single nonzero DFT 
coefficient, while all of its values in the time domain 
are nonzero. That is a kind of signal we say is 
localized in frequency. If a signal is localized in time 
or in frequency, we can efficiently store and process 
its signal values by storing/operating on the signal 
in the best of those two domains. 


But some signals are not localized very well in time 
or in frequency. Supposing a signal has frequency 
components that show up and then are absent at 
different times, then neither the time nor the 
frequency domain does a good job of efficiently 
representing such as signal. An example of such a 


signal would be a musical melody; at any given time 
it may be expressed as being a single (predominant) 
frequency, but these frequencies (i.e., the notes) 
obviously change quite a bit over time (the work of 
Philip Glass notwithstanding). 


There are many transforms that are designed 
specifically for analysis and concise representation 
of these signals that have localized frequency 
components, two of which we will consider. The 
first is the Haar Wavelet, which is a type of wavelet 
transform. The basis functions of the Haar wavelet 
transform, like other wavelets, are localized in both 
time and frequency, as illustrated in the following 
figure: 


; fine scale wavelet 

teeebeteeeb reese 

1 

0 5 10 15 

| mid scale wavelet 

seeeb tpg tloreee. 

1 

0 5 10 15 
coarse scale wavelet 


Unlike the DFT or DCT, wavelet basis functions are 


nonzero for most of their duration in time. But 
unlike delta functions, they have particular 
frequencies for their duration. So in a sense they 
have the benefits of both domains, combined. 

A STFT of a speech sample. 


The Short Time Fourier Transform 


Another transform that considers signals with time- 
localized frequency components is the short-time 
Fourier transform (STFT). The STFT is kind of like a 
DFT, except that it first splits a long signal into 
many smaller parts, and then takes a DFT of each of 
those parts. This "splitting" is done by multiplying 
the signal (of length N) by a moving window (of 
length Nw), progressively taking the DFT of the 
product: 


X[m,k] = Xn =ON — 1x[n]w[n—m]ej2aNwk(n—m). 


The STFT of a signal is therefore usually represented 
as a two-dimensional image. Each column (m) of the 
image corresponds to a time window of a signal, 
while the rows (k) indicate how much of each 
frequency is present at a given time. It is therefore a 
bit like the sheet music score for a piano 
composition, in which the left-to-right axis of the 
music staff corresponds to time, while the markings 
up and down the staff correspond to the notes (i.e., 
frequencies) that are played at the given time. 
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Eigenanalysis of LTI Systems (Infinite-Length 
Signals) 

Complex harmonic sinusoids are eigenvectors of 
discrete-time infinite length LTI systems. When 
given as inputs to any LTI system, the output is 
simply a scaled version of the input. 


Eigenvectors of LTI Systems (infinite- 
length) 


For a system H operating on infinite-length discrete- 
time signals, an eigenvector is a signal which, when 
passing through the system, is modified only by a 
scalar factor. So suppose s[n] is an eigenvector to 
system H, then H{s[n]} =Ass[n], where As is a 
complex number. It is a remarkable property that all 
LTI systems have a set of eigenvectors, and 
furthermore that these eigenvectors are of the form 
ejan. This means that, whatever w may be, if ejan 
is input into any LTI system, the output will simply 
be a scaled version of the input (the amount of 
scaling depending on the nature of the system, of 
course): 
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Proving this special property of LTI systems is 
straightforward; we simply express the output in 
terms of the complex harmonic sinusoid input 
(soo[n] =ejoon) and system impulse response (h[n]) 
through the convolution sum, and simplify: 


so[n]*h[n]=X%mMmM=—-—-x ocswo[n—-m] 
h[m]=2m=- © eejo(n-—-m)h[m] = 
Xm=-ecejone—-jomh[m]=(2m 


=-—-cecoh[mle-—jom)ejan=Awsal[ 
n | 


The Frequency Response of an LTI System 
Take a close look at the value Aw above. It is 
intimately related to the system, as we would 


expect, by the system's impulse response: 


A. = XM = — © oh[mJe—jom. 


This is simply the DTFT of the impulse response! It 
is so important, we give it a special notation H(a), 
and a special name as well. H(@) is called the 
frequency response, for it explains how the system 
responds (i.e., in what way it scales) to particular 
input frequencies w. This is also apparent by seeing 
the DTFT as an inner product of h[n] with the 
sinusoidal e—jam. The value of the inner product 
(which is the value H(o) grows or shrink in ways 
corresponding to the similarity of h[n] and e—jom). 


Although we cannot display the system matrix 
involved (as it is infinite in length), the DTFT 
diagonalizes the system H for infinite length signals, 
just as the DFT does with finite length ones. 
Expressing this diagonalization in terms of the 
DTFTs of the system, we have that if 

yn] =h[n]*x[n], then Y(w) = H(o)X(o). 


The Discrete-Time Fourier Transform 


It is a significant concept in signal processing that 
there is more than one way to express the 
information in a given signal. One way is by 
defining the signal in terms of the value it has at 
any given time n, thus the notation x[n]. But is also 
possible to express a signal in terms of how it is be a 
combination of signals from a particular signal set. 
We could express the information in the signal by 
noting how much of each signal in the signal set 
contributes in the combination. 

The DTFT of an infinite-length signal can be 
thought of as the limit of adding an infinite number 
of zeroes before and after a finite-length signal. 


Definition of the DTFT 


The discrete-time Fourier transform, or DTFT, is a 
common way to express infinite-length discrete-time 
signals in another way. As noted above, the signal 
x[n] can be understood as an entity as having a 
particular value for every time n in the range 

— co <n< oo, But it also can be understood in this 
way: 


x[n]=12nf—-xxX(wo)ejonda 


x[n] is expressed in terms of how much of each 
signal ejwn contributes to the integral which 


composes it. The weighting of each signal is 
according to the function X(@). The values of that 
function are found according to this formula: 


X(@) = Xn = — o ox[nle—jon. 


This formula for X(w) is known as the DTFT of x[n], 
while that first integral formula is known as the 
inverse DTFT. 


Interpretations of the DTFT 


We have given the formulas for the DTFT and the 
inverse DTFT: 


*X(o)=X%n=-e7e ox[nle-jon 
*x[n]=12nf—-axaxX(@a)ejondw 


But what do these formulas mean? There are a few 
ways of understanding them. One is to say that the 
two formulas are expressing the same entity, but in 
different domains: x[n] is the time-domain 
representation of the signal, whereas X(w) is the 
frequency-domain representation of it. In that sense, 
it is like two different language translations of the 
same book (this analogy is not perfect, for 
translations are not exactly one-to-one; there are 
many English-language translations of the one 
Spanish book Don Quixote, but there is only one 
frequency-domain representation X(@) for a time- 
domain signal x[n]). A second way to understand 


the DTFT is in terms of the concepts of analysis and 
synthesis. The DTFT formula gives us a frequency 
analysis of the signal x[n]. The signal X(@) analyzes 
x[n] from a different point of view; it tells us the 
frequency content of x[n], whether it is composed of 
only low frequencies, or just high ones, or whatever 
the case may be. This is because every value of X(o) 
is simply the inner product of x[n] with ejwn; it tells 
us the "strength" of ejon in the signal x[n]. The 
inverse DTFT then synthesizes, or composes, x[n] 
based upon its weight at each different frequency. 


Relationship between the DTFT and the DFT 


The DTFT is not the only Fourier transform for 
discrete-time signals. While the DTFT operates on 
infinite-length signals, the DFT (discrete Fourier 
transform) is a tool for frequency analysis of finite- 
length signals. A finite-length signal (say, of length 
N) x[n] has an N length frequency representation 
X[k], according to these formulas: 


X[k] = xn = ON — 1x[nJle —j2aNknx[n] = 1NXk=0ON 
—1X[k]ej2xNkn. 


As x[n] and X[k] are N-periodic, we can also 
represent them in this way, shifted about n and k at 
the origin: 


X[k] = xn = —N2N2 — 1x[nJe—j2aNknx[n] = 1INZk= 
—N2N2—1X[k]ej2xNkn. 


We could develop the DTFT by seeing what happens 
to the DFT as N gets larger and larger. As N 
approaches -, so also do n and k run from — © to 
co, Consider also what happens to e—j2aNkn. As N 
gets infinitely large, the exponent term 2Nk 
becomes a continuous value, running from — zx to 7, 
and we can call that value ow. We have then that the 
DFT sum becomes: X(@) = Xn = —  ox[nle—jon 
which is the definition of the DTFT. Similarly, on 
the synthesis side of the transform, the DFT "sum" 
will not run over a discrete number of values k, but 
rather become an integral over the continuous 
variable o: x[n] = { —1X(w)ejondw2x which is our 
inverse DTFT. Graphically, we can see how 
increasing the length of a signal by adding more and 
more zeros on either side of it results in the signal's 
DFT (which is expressed in terms of k) appearing to 
be more and more like a continuous function 
(expressed in terms of w): 
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Discrete-Time Fourier Transform Examples 


The DTFT and inverse DTFT are defined as follows: 


X(w) = Un=—-eox[nle—-jon, 
—-nsom<aux[n] = J[-—nuxX(@)ejo 
ndw2za, coo <n < © 


Let's work out some examples of DTFT and inverse 
DTFT calculations. 


Impulse Response of an Ideal Lowpass 
Filter 


We'll start with an ideal lowpass filter. From its 
frequency response we can see that it blocks all 
incoming frequencies having a magnitude greater 
than |ac|: 


H(o) = {1 —wce <= < wcOotherwise. 


In order to find the impulse response for such a 
filter, we will perform an inverse DTFT on the 
frequency response: 


h[n]=jfJ—-xxnH(o)ejondw2n=122f 
-ococlejondwo =12x[ljnejon] — 
ococ=l1rn(ejwocn—ej@cn)2j =sin ( 
@cn)mn 


This impulse response is a form of the sinc function, 
and is plotted below. 
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DTFT of a Moving Average System 


Consider now a moving average system, where the 

system output at a time n is the average of M input 
values, symmetric about n. The impulse response of 
such a system is a pulse function: 


p[n] = {1 -M<n<MOotherwise. 
p(n 


Suppose we would like to find the frequency 
response of this system. We merely need to take the 


DTFT of the impulse response: 


P(wm) = Xn= — ~ op[nje—jon = Xn= —MMe 
—jon = Yn= —MM(e—j@)n = ejaM—e—jo(M 
+1)1—e—-jw = ejmM—e—-jo(M+1)1—e—-jwo =e 
—jo/2(eja2M + 12—e—jw2M+12)e 
—jo/2(eja/2—e—j@/2) = 2jsin(w2M 
+12)2jsin(w2) = sin(2M+12)sin(w2). 


This frequency-domain representation of the time- 
domain pulse is known as a digital sinc: 
P(w) 


The digital sinc function is similar to a regular sinc 
function, except that rather than perpetually 
decaying, it is 2% periodic. When we compare this to 
the first example, the impulse response of an ideal 
lowpass filter, we see an interesting correspondence. 
Pulses in one domain (either time or frequency) 
have sinc-like representations in the alternate 
domain. 

The impulse response of a recursive average 
system. The magnitude of the frequency response of 
a one-sided exponential. 


DTFT of a One-Sided Exponential 


Having found the frequency response of a moving 
average system, let's now do the same for a 
recursive average system. The input/output 
relationship of such a system can be expressed as 
y[n] =x[n] + ay[n—1], where |a| <1. The impulse 
response for such a system is h[n] =anu[n]: 

h[n] = a" u[n], a=0.8 


We'll now take the DTFT of this impulse response: 


H(o) = Xn= — © oh[nje—jon = Xn=O-~ane 
—jon = Yn=0~(ae—jo)n = 11—ae—jo. 
|H(w)| 


The Discrete-Time Fourier Transform of a Sinusoid 
Recall the definition of the DTFT of a signal x[n]: 
X(@) = XYn= — ~ ox[nje—jon. 


Let's see what the DTFT of a complex exponential 
signal ejwOn: 


X(w0) = Xn= — ~ ox[nle—jon = Yn= 
— co cej@One—jon = Xn= — © e“e—j(@—w0)n. 


The nature of that sum is not immediately clear. 
When w=0, then the sum becomes unbounded: 


Yn=—-wol=ao, 


But when w ~#w0... what does an infinite sum for a 
sinusoid mean? 


The Dirac Delta Function 


To find the DTFT of a sinusoid, it will to be 
necessary to consider another function, the dirac 
delta function. Suppose we have a function de(): 


Note how the non-zero portion of this function has a 
width of e and a height of 1e. The area underneath 
the function is therefore always 1, no matter what € 
may be. We will let the value of € get smaller and 
smaller. As it approaches zero, the function will 
become taller and narrower, until at the limiting 
case it is "infinitely" tall and "infinitesimally" 
narrow. We will denote the function in this limiting 
case 6(). But it will still have an area of 1, thus 


fd(o)dw=1. 


Suppose we scale the 6(@) by another function 
X(o). Since 5(@) is non-zero at all values other than 
at w=0, scaling 6(@) by X(@) within the integral 
will yield the following: 


JX(w)5(@)dw = X(0). 
What if the integrand has 5(w — 0) and is scaled by 
X(o). In this case the integrand is only nonzero 


at wo = 0, thus: 


JX(w)5(@ — 00)da = X(w0). 


The 6(o) "function" we have derived is known as 
the Dirac delta function. It is not technically a 
function (it is a generalized function, or 
distribution), but for our purposes we can treat it as 
a function with infinite height, infinitesimal 
narrowness, and an area of 1 at w=0. 


The DTFT of a Sinusoid 


Having introduced the Dirac delta function, let's see 
what happens if we attempt to find the inverse 
DTFT of 23t8(@ — w0): 


J —11216(0 — oO)ejandw2n = ejwOn. 

The inverse DTFT of a dirac delta function is a 
complex sinusoid, which means that the DTFT of a 
complex sinusoid is a dirac delta: 


eja0n <-DTFT 275(w — w0). 


Via Euler's formula, we use the above to also find 
the DTFT of discrete-time cosines and sines: 


cos(wOn) <-DTFT 18(@ — 00) + 15(@ 
+ w0)sin(o0n) <-DTFT rj5(@ — w0) + tj5(w + w0). 


Discrete-Time Fourier Transform Properties 


Recall the definition of the DTFT and the inverse 
DTFT: 


X(m) = Xn= — co ox[nJe 
—jon, —t=w<ax[n] = f 
—mX(w)ejondw2r, coo <N<oo , 


A signal and its corresponding DTFT is referred to as 
a DTFT signal pair: 


x[n] ©DTFT X(o) 


We'll now take a look at some properties of the 
DTFT. 
The 2x periodicity of the DTFT. 


DTFT Periodicity 


Although the DTFT of a signal is typically defined 
only on a 2x interval, it is nevertheless 
straightforward to show that the DTFT is actually 27 
periodic: 


X(@+ 20k) = Xn= — ~ ox[nle—j(w 
+ 2nk)n = Xn= — © ox[nle—jwne—j2rkn = X(o). 


DTFT Frequencies 


Because X(w) is essentially the inner product of a 
signal x[n] with the signal eja@n, we can say that 
X(@) tells us how strongly the signal ejon appears 
in x[n]. X(w), then, is a measure of the "frequency 
content" of the signal x[n]. Consider the plot below 


of the DTFT of some signal x[n]: 
X (w) 


This plot shows us that the signal x[n] has a 
significant amount of low-frequency content 
(frequencies around w= 0), and less high-frequency 
content (frequencies around w= +12 -- remember 
that the DTFT is 22 periodic). 


As the DTFT is 2x periodic, all the information in a 
DTFT is contained in any 2s length section. 
However, there are only two intervals that are 
typically used when representing a DFTF of a signal. 


The first is to use the interval of —2 to , as above. 
In such an interval, the low frequencies are at the 
center, with the high frequencies at the edges. The 
other interval ranges from 0 to 2. The low 
frequencies are at the extremities of the plot, with 


the highest frequency in the center: 
X(w) 


The DTFT and Time Shifts 


If a signal is shifted in time, what effect might this 
have on its DIFT? Supposing x[n] and X(o) area 
DTFT pair, we have that 


x[n-—m] <DTFT e—jamX(@). 


So shifting a signal in time corresponds to a 
modulation (multiplication by a complex sinusoid) 
in frequency. We can use the DTFT formula to prove 
this relationship, by way of a change of variables 
r=n-—m: 


Mn=-eeox[n-mjle-jon=2r= —- 
co ox[rle-j@a(r+m)=Yr=—-eo~x|[ 
rle-jore-—jom=e-jom2Zr= —- © » 


x[rle-jwr=e-—-jomx(o) 


The DTFT and Time Modulation 


We saw above how a shift in time corresponds to 
modulation in frequency. What do you suppose 
happens when a signal is modulated in time? If you 
guessed that it is shifted in frequency, you're right! 
If a signal x[n] has a DTFT of X(o), then we have 
this DTFT pair: 


ejwOnx[n] <?DTFT X(@— 0). 
Below is the proof: 


rn=e cejoOnx[n]e-—jwnun= —- ~ 
ox[nle-j(a—-oa0)n=X(ow-o0) 


The DTFT and Convolution 


Suppose that the impulse response of an LTI system 
is h[n], the input to the system is x[n], and the 
output is y[n]. Because the system is LTI, these 
three signals have a special relationship: 


y[n] = x[n]*h[n] = Xm= — ~ oh[n—m]x[m]. 


The output y[n] is the convolution of x[n] with 
h[n]. Just as with the other DTFT properties, it 


turns out there is also a relationship in the 
frequency domain. Consider the DTFT of each of 
those signals; call them H(@), X(@), and Y(@). The 
convolution of the signals x and h in time 
corresponds to the multiplication of their DTFTs in 
frequency: 


y[n] =x[n]*h[n]«?DTFTY(@) = X(@)H(). 


For the proof, we take the DTFT of y[n], using a 
change of variables along the way: 


Y(o)=%n=- er wey[n] =AnN= — ~ ~ ( 
x[n]*#h[n])e-jon=XZn= —-e2 om 

=-—-eoox[m]h[n-—-m]le-jon=%Yym= 
-ooYn=—-—-eooxl[m]lh[n-—-mle-jo 
me-jo(n-—-m)=XmM=-% ox[mle- 
jom(rn=—-e oh[n-—-mjle-jo(n-m 
)) Let V=en-—-m=2X%M=—-er%ox[mle-j 
am(Xv=-*7 -meo —-mh[vle-jo(v) 
J=Xim=—-er ex[mle-—jom(Xuv= —- ~ 
ch[vle-—jw(v))=Xm= —-eox[m]le 
—-jom(H(o))=H(o)xem= —-©7 ox[m 


le-joam=xX(o)H(o) 


This relationship is very important. It gives insight, 
showing us how LTI systems modify the frequencies 
of input signals. It is also useful, because it gives us 
an alternative way of finding the output of a system. 
We could take the DTFTs of the input and impulse 
response, multiply them together, and then take the 


inverse DTFT of the result to find the output. There 
are some cases where this process might be easier 
than finding the convolution sum. 


The DTFT and Linearity 


Since the DTFT is an infinite sum, it should come as 
no surprise that it is a linear operator. Nevertheless, 
it is a helpful property to know. Suppose we have 
the following DTFT pairs: 


xl[n] <-DTFT X1(o@) x2[n] <?DTFT X2(o). 


Then by the linearity of the DTFT we have that, for 
any constants a1 and a2: 


alxifn]+a2x[2] -—DTFT alX1(o 
y+t+a2X2(o) 


DTFT Symmetry 


Looking at the plots of the DTFTs above, or perhaps 
of those you have computed on your own, you may 
notice that they are often symmetrical. It so happens 
that the DTFTs of real-valued signals always have a 
kind of symmetry, but there are other symmetry 
relationships as well. For example, signals that are 
purely imaginary and odd have DTFTs that are 


purely real and odd. These types of symmetry are a 
result of a property of the complex exponentials 
which build up the DTFTs. Any signal of the 

form ej@n is conjugate symmetric, meaning that its 
real part is even and its imaginary part is odd. 
Additionally, for conjugate symmetric signals, their 
magnitude is even and their phase is odd. 


As a result of the symmetry of the signals which 
compose a DTFT, the DTFTs of certain signals also 
have symmetries, according to the table below: 
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The z-Transform 


Generalizing the DTFT 


We have seen that the DTFT is a useful tool for 
representing and analyzing discrete-time infinite 
length signals and LTI systems. With the DTFT, 
signals can be represented in terms of their 
frequency composition. Because complex harmonic 
sinusoids are eigenvectors of LTI systems, the DTFT 
of a system's impulse response gives us the system's 
frequency response (it tells us how the system 
modifies input signals according to their composite 
frequencies). 


The backbone of the DTFT is the complex harmonic 
sinusoid: ejon. This type of signal is actually a 
specific example of a more broader class of signals 
zn, where z is a complex number. Signals in this 
class are called complex exponentials. As zn = |z| 
nej(2z)n, we see that complex exponentials are 
simply complex harmonic sinusoids, with 
exponential function envelopes. When |z| <1, they 
decay in time: 
Im(z"), |z| <1 
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And when |z| =1, they are of course simply complex 
harmonic sinusoids of the form ejon. 


Since complex harmonic sinusoids have the special 
property that they are eigenvectors of LTI systems, 
we may ask if the more general signal class of 
complex exponentials does, as well. To test, we will 
find the an LTI system's output (i.e., the result of the 
convolution sum) to some generic complex 
exponential zn: 


yn] =zn*h[n] = Xm= — ~ ox[n—mJh[n] =xym= 
— co cooz7n—mh[n] =xym= — ~ ©Znz 
—mh[n] =znim= — © ©z—mh[n] =Qm= — © 07 


—mh[n])zn=H(z)zn , H(z) = Xm= — © ©z—mh[n]. 


If we input a complex exponential into some LTI 
system, then indeed the output is simply a scaled 
version of that same exponential. So complex 
exponentials are also eigenvectors of LTI systems. 
Now, the amount a given LTI system scales a 
complex exponential zn is very important, so much 
so that it has a special name--it is the system's 
transfer function (it tells us how the system 
"transfers" the input into the output)--and a special 
label, H(z). 


This transfer function of LTI systems, H(z), is clearly 
a generalization of the frequency response H(o). 
The frequency response indicated how an LTI 
system H scales input complex sinusoids ejw, and 
likewise the transfer function H(z) tells us how the 
system scales input complex exponentials zn. As 
ejon = zn|z = ejw it follows that H(@) = H(z)|z=ejo. 
So it would seem that the DTFT of a function (for 
H(@) is the DTFT of h[n]) can be generalized to be a 
broader kind of transform, which is represented by 
H(z). This broader class of transform is called the z- 
transform. 


The z-Transform 


Let x[n] be a discrete-time infinite length signal. 
Then X(z) = Xn= — ~ x[n]z—n is defined to be the 


bilateral z-transform of x[n] (the unilateral z- 
transform sums only from n=0 to ~; unless 
otherwise stated, we assume the transform is 
bilateral). 


There are a few things we can note immediately 
about the z-transform. First, unlike either transform 
we have considered (the DFT and DTFT), it is a 
complex function of a complex variable. Second, and 
like the DTFT (but not the DFT), the z-transform 
may not necessarily exist for all z and/or all x[n]. 
And third, like both other transforms, the z- 
transform diagonalizes LTI systems. By this, we 
mean that if a signal x[n] is expressed in terms of its 
z-transform, then the z-transform of the output at 
various z is simply the pointwise multiplication of 
the input's transform with the transform of the 
impulse response (i.e., the transfer function: 

Y(z) = H(z)X(z)). 


The z-Transform Region of Convergence 


The Existence of the z-Transform 


Recall the definition of the z-transform for some 
discrete-time infinite length signal x[n]: 


X(z) = Xn = — c ox[n]z—n. 


Given that this is an infinite sum (unlike, say, the 
sum of N terms in a DFT), it will not necessarily 
converge for all x[n] and/or all z. For some x[n], it 
will converge for all values of z. For example, 

if x[n] =5[n], then: 


X(z) =Xn= — © o§[nJz—n=1. 


For some x[n], the sum will not converge for any z. 
For example, if x[n] =1, there is no z such that the 
sum converges. 


The Region of Convergence 


So for some x[n] (like x[n] =6[n]) the z-transform 
sum will converge for all z, and for some x[n] (like 
x[n] = 1) the z-transform will not converge for any z. 
As you might suspect, there are some x[n] for which 
the sum converges for some z. These values of z, 


those for which the z-transform exists, are known as 
the z-transform's region of convergence. Consider, 
for example the z-transform for the signal 
xl[n]=anu[n] , a=.8: 

x(n] =a" u[n], a=0.8 


X1(z) =Xn= — © ox1[n]z—n=Yn=0-~-anz 
—n=Xn=0~(az—1)n={1l—az—-1laz-1| 
< 1undefinedelse. 


So, the z-transform X1(z) only exists for |az—1| <1, 
or put another way, for |z| >a; it is not defined for 
other values of z. This should not be too terribly 
shocking, for we are familiar with other functions 
that only exist for certain values, e.g. logx only for 
x >0, or 1x only for x +0. 


However, there is a difference between the z- 
transform and its region of convergence and 
functions like logx or 1x. Where the defined regions 
of existence of the values of logx or 1x are quite 
clear, there is nothing intrinsic about the function zz 
—a to suggest it only exists for |z| >a (of course it 
is apparent that z cannot equal a). Why in the world 
would it be a problem if z=a/2, for then we would 
simply have a/2a/2—a= —1, right? The reason is 
that the function was the simplification of the sum, a 


simplification that existed only for certain z. 


To see why that is the case, consider another signal 
related to the one we just considered. Let x2[n] = 
—anu[ —n—- cr a=.8: 

ro[n] = —a~" u[—n — I}, a=0.8 
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We will now consider its z-transform: 


X2(z)=%n=—-—-ecox2[n]z-—-n=Xn= 
= = 1S ON nS SS bea 
=-'Yn=Oe(a-—-lz)n+1={-ll-a 
—1z+1]a-—1z| < 1 undefined else = { -— 1 
DS 2h a 2 ee 2 ord 
z| < 1 undefined else = {11—az-—lla-—1l1z 
| < 1 undefined else 


So for x2[n], its z-transform is defined only for |a 
—1z| <1, or equivalently, |z| <a, and for those 
values it is 11-—az—1. Now we see another reason 
why the region of convergence is so important; it is 
not merely added information about the z- 
transform, it is intrinsic to its very definition. The z- 
transforms for x1[n] and x2[n] are completely 
different, which is to be expected as those are 
completely different signals. However, the only 
difference between the transforms is the region of 


convergence: 


xl[n]=anu[n]?~X1(z)={1ll—-—az- 
1|z| > a undefined elsex2[n] = -—anu[- 
n-1J<+X2(z)={ll-—az-1l1|z|<a 

undefined else 


Types of ROCs 


We have seen a few different types of ROCs for the 
z-transforms of various signals. If the z-transform for 
a discrete-time signal exists, then the region of 
convergence will fall into one of these three 
categories: 


¢ If the signal has a finite non-zero duration, then 
the region of convergence will be the entire z- 
plane, with the possible exception of z= © (if 
the signal has nonzero values at time indices 
less than zero) or z=0 (if the signal has 
nonzero values at time indices greater than 
Zero). 

¢ If the signal is of infinite duration as n 
approaches - (but not — ~), then the region 
of convergence will extend outside some disk 
(e.g., |z| >1) 

¢ If the signal is of infinite duration as n 
approaches — - (but not ~), then the region 
of convergence will be inside some disk (e.g., | 
Z| <4) 


¢ If the signal is of infinite duration as n 
approaches both ~ and — ~, then the region 
of convergence--if it exists at all--will be an 
annulus (e.g., 1 < |z| <2) 


The Transfer Function of Discrete-Time LTI Systems 


LTI Systems in the Time Domain 


One of the most general ways of defining an LTI 
system is by expressing its output at some time n as 
a weighted sum of values of the input and output at 
a finite number various other times. If all those 
times are at or prior to n, such a system is causal. In 
terms of a difference equation, the system definition 
would look something like this: 


y[n] = bOx[n] +blx[n-1]+b2x[ 
n-2]+-°:+bMx[n-—-M]-aly[n-1] 
a? Ya es Do) eo = aN PN 


This kind of system can also be displayed as a block 
diagram: 


x[7] 


In the block diagram, the z—1 blocks represent time 
delays, the triangles represent multiplications, and 
the & blocks represent summations of the signals 
being input into them. The output at time n is a 
weighted sum of the previous M input values (along 
with the current value) and previous N output 
values. 


LTI Systems in the z-Domain 


Consider again the input/output relationship of a 
causal LTI system: 


y(n] = bOx[n] +b1x[n—1]+b2x[n-—2]+-::- 
+ bMx[n—M] —aly[n-1] —a2y[n—2] +-+-—aNy[n 
—N]. 


We can take the z-transform of this whole equation, 
bearing in mind that the transform is a linear 
operation and that delays of nO in time correspond 
to multiplications by z—n0 in the z-domain. That 
would leave us with this: 


Y(z)= bOX(z) + blz—1X(z) + b2z — 2X(z) + ---bMz 
—MX(z)-—alz—1Y(z) —a2z—2Y(z) —-:-—aNz 
—NY(z). 


From there we can do a little rearranging: 


Y(z) = bOX(z) + b1z — 1X(z) + b2z — 2X(z) + -:-bMz 


—MX(z)-—alz—1Y(z) —a2z—2Y(z) —---—aNz 
—NY(z)Y(z) + alz—1Y(z)a2z — 2Y(z) +--+: +aNz 
—NY(z) = bOX(z) + b1z—1X(z) + b2z— 2X(z) + --bMz 
—MX(z)¥(z)(. +alz—1+a2z—2+---+aNz 
—N)=X(z)(b0 + blz—1+b2z—2+---bMz 

—M)H(z) = Y(z)X(z) = b0 + blz—1+b2z—2+---bMz 
—Ml1+alz—1+a2z—-2+::-+aNz—N. 


This H(z), which is a rational function in terms of z, 
is called the system's transfer function. 

A system with two poles has three possible regions 
of convergence.Consider a system with two poles, as 
plotted here. One region of convergence could be a 
disc within the smallest magnitude pole. Another 
region of convergence could be an annulus between 
the two poles. A final region of convergence could 
be all values of z of greater magnitude than that of 
the largest pole. 


The Transfer Function: Poles and Zeroes, 
ROG, and Stability 


So we see that the transfer function for one of the 
most common classes of LTI systems--in which the 
output is a weighted sum of a finite number of past 
inputs and outputs--is a rational function in terms of 
vs 


H(z) =b0+bl1z—1+b2z-—2+-:-bMz—M1+alz 
—1+a2z7-—2+°+aNz—N. 


We can express this function in terms of positive 
powers of z by simple factoring, and then according 
to the fundamental theorem of algebra, we can 
represent the numerator and denominator 
polynomials with respect to their roots: 


H(z)=b0+bl1z—1+b2z-—2+::-bMz—Ml1+alz 
—1+a2z-—2+:-:-+aNz—N=z—M(b0zM+ b1zM 
—1+b2zM—2+::-+bM)z—N(zN+a1zN —1+a2zN 
—2+:--+aN)=z—Mb0(z— (1)(z— €2):--(z-— CM) z 
—N(z— p1)(@— p2)---(z— pN). 


The values of z for which the numerator is zero are 
called zeros, and the values for which the 
denominator is zero are called poles. 


Now, there is a special reason the transfer function 
is given the notation H(z): it happens to be the z- 
transform of h[n], the impulse response of the 
system. Like all z-transforms, it has a region of 
convergence. We have seen previously that the ROC 
for a z-transform will be one of the following: 


* no z (the z-transform does not converge for any 
Z) 

* the entire z plane (with the possible exception 
of z=0 or z=) 

- all z outside a certain disc (e.g., |z| >1) 

- all z within a certain disc (e.g., |z| <1) 

* z within an annulus (e.g., 1<|z| <2) 


Poles play a very significant role in the ROC for a 


transfer function; in fact, they completely define the 
ROC. The ROC of a z-transform cannot contain a 
pole (for the function would not exist at that point), 
which means that once we know the poles of a 
system, we know the potential ROCs. We say 
potential because there are a variety of possible 
ROCs for a transfer function with one or more poles; 
to be precise, if there are M distinct poles, then 
there will be M+1 possible ROCs: the disc within 
the pole of smallest magnitude, the region outside 
the disc or radius equal to the largest pole, and then 
the M—1 annuli of regions between adjacent poles. 
If, suppose, the system in question has two poles, 
then there will be three possible regions of 
convergence: 
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As we have also seen previously, each type of ROC 
corresponds to a particular kind of discrete-time 
signal. ROCs which are a disc correspond to signals 
with a non-zero duration extending to n= — ~, 
those that are all z outside of a certain disc 
correspond to signals with a non-zero duration 
extending to n= ~, and those that are an annulus 
correspond to signals whose non-zero duration 
extends both to n= — ~ and n=. For the class of 
discrete-time LTI system that we have been just 
considering, the impulse response is causal, and thus 
if it is infinite, will have a duration that extends 
towards n= ©. Thus the transfer function of such a 
system has a ROC that extends outward from the 
outermost pole. 


Finding the Transfer Funcion, Zeros and Poles, 
and ROC 


Consider a discrete-time system represented by this 
block diagram: 


x(n] [2 Dales are 


ne a\ EN aN 


=| ay 


We would like to find its transfer function, along 
with its poles and zeros. The first step will be to 
express the input/output relationship in the time 
domain. From there we can take the z-transform of 
each side and find Y(z)X(z), which is by definition 
the transfer function: 

y(n] =x[n] + 3x[n—-1]+114x[n—2]+34x[n 

Soli yt et yay a ya et yi 
—2|=x([n] +3x[n—1]+114x[n—2]+34x[n 

— 3]Y(z) —z—-—1Y(z) +12z —2Y(z) = X(z) + 3z 
—1X(z)+114z—2X(z) + 34z-—3X(z)Y(z)Q —-z 

Se 2) le or i ae 

— 3)H(z) = Y(z)X(z) =14+3z-14+114z-24+34z 
—31-—z-14+12z-—2=2-3(23+3z24+114z+4+ 34)z 
— 2(z22 —z2+12)=(z24+1)(z+12)(Z2+32)z7 
=A Ep) 

So, having found the transfer function and factored 
the numerator and denominator, we see readily 


that the zeros of the system are z= —1/2, —1, 

— 3/2, and the poles for the system are z=0, 
1/2+j/2, 1/2 —j/2. As the system is causal (from 
the block diagram we see the output is dependent 
only on present and/or past values of the input and 
output), the ROC will extend outward from the 
outermost pole. The two nonzero poles are 
equidistant from the origin, at a distance of 22, so 
the ROC is |z| >22. A pole/zero plot, with the ROC 
shaded, is below: 


Im(z) 


From the Transfer Function to the 


Frequency Response 


We have seen that for LTI systems whose output is a 
weighted sum of a finite number of past inputs and 
outputs, the transfer function amounts to being a 
rational function of the form H(z) =b0 +b1z 

=1 b2z]24"DMz= M1 4+alz—1a2z—24 
+aNz—N, or equivalently, z— Mb0(z— C1)(z—¢2)::: 
(z—CM)z— N(@-— p1)(z— p2)::-(z— pN). One of the 
benefits in expressing the transfer function as the 
products of those single order polynomials is the 
elegant way it leads to understanding the magnitude 
of the transfer function: 


|H(z)| = |z—M||bO||(z—¢1)||(z—¢2)|---|(@— CM) ||z 
—N||(z—p1)||(@@—p2)|---|(z2—pN)|. 


Note what terms like |z—p| represent. The 
magnitude is (Re (z) — Re (p))2+ Um (z) —Im (p))2, 
which is simply the euclidean distance between z 
and p on the complex plane. So the magnitude of 
the transfer function at some point z essentially 
amounts to being the products of the distances 
between z and the system's zeros, divided by the 
products of the distances between z and the system's 
poles. The closer z is to a zero, the smaller the 
magnitude of the transfer function will be; the closer 
z is to a pole, the larger the magnitude of the 
transfer function will be. 


Now recall that a system's frequency response, i.e., 


the amount the system scales the input at particular 
frequencies, is simply the z-transform evaluated on 
the unit circle, for H(ejm) = H(z)|z=ejw. If we have 
a pole-zero plot, we can quickly visualize the 
magnitude of the system's frequency response by 
traversing around the unit circle. As we go, the 
closer we are to a zero, the smaller it will be, while 
the closer we are to a pole, the larger it will be. As a 
result of this phenomenon (and the fact that the 
transfer function is determined by the pole and zero 
locations), we can see that the locations of a 
system's poles and zeros in the z-plane totally 
characterize a system (within a scaling factor). Thus, 
knowing the poles and zeros of a system is 
essentially equivalent to knowing exactly what kind 
of system it is. 


isualizing the Magnitude of the Frequency 
Response 
Consider the pole-zero plot from the example 
above, with a unit circle overlaid: 


Im(z) 


Re(z) 
The frequency response of this system is simply 
value of the transfer function, evaluated along this 
unit circle. The magnitude of the transfer function 
is the product of the distances to the zeros, divided 
by the product of the distances to the poles. By 
traversing along the unit circle and noting these 
distances, we can visualize what the magnitude of 
the frequency response will be. We will start at 
«= 0. Here we can easily find the value of the 
magnitude, since H(ej0) =H(1). From the transfer 
function we found in the first example, we have: 
H(1)=(1+34+114+34)0 -—14+12)=15. 
So |H(ej@)| at o=0 is 15. Now, as we move 
counter-clockwise around the circle (i.e., from 
w=0 to o=12 to ®=7), note what happens. As w 
approaches 214 on the unit circle, it is drawing 
nearer and nearer to a pole; therefore, the 
magnitude of the frequency response is going to 


increase as w progresses from 0 to 114. At that 
point, we will be travelling away from that pole, 
and nearer to the three zeros; accordingly, the 
magnitude of the frequency response will decrease. 
Eventually we will actually land squarely on a zero 
at «=m; this is obviously as close as we could 
possibly get to a zero. As the distance to the zero at 
that point is zero, the magnitude of the frequency 
response there is zero. We can repeat this process 
moving clockwise from w@=0 to w= —z and the 
magnitude will change in the precisely same way, 
due to the symmetry. Having moved along the 
entire unit circle, we can see how our intuition 
matches what is indeed the magnitude of the 
frequency response, found via plotting software: 


Note that the magnitude is 15 at w=0, then 
increases to a peak at about +714, and then 
decreases to 0 as w traverses from +74 to +n. 


The Transfer Function ROC, and Stability 


It has been explained previously that the region of 
convergence of a z-transform are the z for which the 
sum »x[n]z—n converges. There is an alternative, 
though related definition that is also often used in 
signal processing, that the region of convergence are 
z for which the sum %x[n]z—n converges absolutely, 
i.e. that &|x[n]z—n]. For virtually all of the signals 
that we will be considering, those two definitions 
are equivalent. The reason the second definition is 
often used is because of a relationship it establishes 
with the stability of the system. Suppose that a 
system's transfer function H(z) has a particular ROC 
(defined in the absolutely summable sense). If the 
unit circle (|z|=1) is contained by this ROC, then 
the system is BIBO stable. To prove this fact, we see 
simply that the |z|=1 being in the ROC means that 
~|h[n]ejw| converges for all w, and then note that 
for @=0, the convergent sum is 2|h[n]|. Previously 
we proved that for LTI systems, %|h[n]] converging 
is equivalent to the system being BIBO stable. 


So if the unit circle is in a transfer function's ROC, 
then that system is BIBO stable. This is a helpful 
way to determine stability. We do not need to prove 
that an arbitrary bounded input will produce a 
bounded output, nor find the impulse response and 
determine its absolute summability; we simply need 
to know the transfer function's ROC. If the system in 
question is causal (which is the case for most all 
practical systems of interest), then the ROC will 
extend outwards from the outermost pole. This 


means that, for causal systems, the location of the 
poles also determine BIBO stability. If all of a causal 
system's poles lie within the unit circle, then the 
ROC will extend to include the unit circle, and 
hence the system will be BIBO stable. 


Stability and the Pole/Zero Plot 
Consider again the pole/zero plot from the system 
in the previous examples: 


Rel z) 
Since the ROC contains the |z|=1 unit circle, this 
system is BIBO stable. 


Poles/Zeros and FIR/IIR Systems 


We have seen that the poles and zeros of a system 
characterize it, which of course includes the nature 
of its frequency response and whether or not the 
system is BIBO stable. Given how significant the 
poles and zeros of a system are, it should come as 
no surprise that knowing about the poles and zeros 
also allows us to determine whether the system is 
FIR or IIR. Put simply, if a system's transfer function 
has any poles (other than at z=0), then the system 
is IIR. If it has only zeros, then it is FIR. This follows 
from what we saw with regard to the nature of 
different kinds of ROCs. If the ROC is the entire z 
plane (with the possible exception of z=0), then the 
impulse response is of finite duration; hence the 
system is FIR. But if there are any non-zero poles on 
the plane, than the ROC will either be a disk, an 
annulus, or will extend outside of a disk; in each of 
those cases, the impulse response is of infinite 
duration. 


z-Transform Properties 


As with other signal transforms, the z-transform has 
a number of significant properties, including ways 
in which changes to a signal in one domain impacts 
its representation in another domain. To examine 
these properties, we will again use the notation of a 
z-transform "pair." If we have X(z) =xn= 

— co cox[n]z—n, then we express the relationship 
between x[n] and its z-transform X(z) as x[n]<X(z). 


The z-Transform and the DTFT 


The DTFT is a special case of the z-transform. If we 
have x[n]<>X(z) and the ROC of X(z) contains the 
unit circle, then the DTFT of x[n], X(ejw), is simply 
the z-transform X(z) evaluated on the unit circle 

Z= ejw: 


X(ejw)=X(z)|z=ejo 


The z-Transform is Linear 


Given that the z-transform is simply an infinite sum, 
if follows then that it will be a linear operator. If we 
have x1[n]<X1(z) and x2[n]X2(z), and al,a2€C, 
then: 


alxil[n]+ta2x2[n]<a1l1X1(z)+a2 
X27) 


However, the ROC of the new transform is the 
intersection of those of X1(z) and X2(z): 


ROCal1X1 + a2X2 = ROCX1/)ROCX2. 


The z-Transform and Shifts in Time 


If x[n] and X(z) are a z-transform pair, then: x[n 
—m]<z—mX(z) and the ROC remains the same as 
for X(z) (with the possible addition/removal of a 
pole at zero and/or zero at infinity): 


yn= — © ox[n—m]z—n=2Xr= — ©» ox[r]z—(r 
+m)=2r= — © ox[r]z—rz—-m=z—myr= 
— co cox[r]z—r=z—mxX(z). 


The z-Transform and Modulation 


If x[n] and X(z) are a z-transform pair, then: 
zOnx[n]<>X(zz0). To determine the ROC of the new 
z-transform, substitute zz0 for z in the original ROC 
definition, and simplify. For example, if the ROC of 
X(z) is |z| >1, then the new ROC will be |zz0| >1—| 
z| > |zO|: 


yn = — ce co (ZOnx[n])z-—n=xXn= — © ox[n](Z/ 
zO) —n = X(zz0). 


The z-Transform and Complex 
Conjugation 


If x[n] and X(z) are a z-transform pair, then: 
x*[n]<X*(z*): 


yn= — © oox*[n]z—n=(XnN= — ~ ox[n] 
(z*) —n)* =X*(z*). 


The z-Transform and Time Reversal 


If x[n] and X(z) are a z-transform pair, then: 

x[ —n]<>X(z—1). With the ROC inverted: to find the 
ROG, substitute 1z for z in the original ROC 
expression and simplify: 


n= — co cox[—n]z—n=XyMm= 


— co cox[m]zm= Ym = — © ox[m](z—1)-m=Xfz 
—1). 


The z-Transform and Convolution 


If x[n]<>X(z), h[n]<>H(z), and y[n]<Y(z) all all z- 
transform pairs, and y[n] =x[n]*h[n], then: 


Y(z) =X(z)H(z), with the ROC of Y(z) being the 
intersection of the ROCs of X(z) and H(z). 


Common z-Transform Pairs 


In practice, one does not typically find the z- 
transform for a signal through the calculation of an 
infinite sum, but rather refers to a table to find a 
transform (along with the use of the properties 
listed above). Below is a table of some of the most 
common z-transform pairs: 


x[n]X(z)ROCS[n]llall zu[n]11—-—z- 
1 2) Stale aa) a2 1 2S 1 on 
u[nJll-—az-—1|z|>|/a|-anu[-n- 
1J1l1l-—az-1]|z|<|a| 


The Inverse z-Transform 


The z-transform for a signal x[n] has been defined 
as: 


X(z) = Xn = — co co x[n]z—n. 


As with other transforms such as the DFT and DTFT, 
we would like to recover x[n] from its transform. 
The technical definition for this inverse transform is: 


x[n] = ¢CX(z)zndzj2nz. 


Evaluation of this contour integral requires 
knowledge of complex analysis. If you are handy 
with contour integrals, feel free to evaluate inverse 
transforms directly! But if you aren't, there is 
another way to find the inverse transform, which we 
will now consider. 


Inverse Transform by Inspection 


In order to find the inverse z-transform, we do not 
necessarily need to calculate anything; rather, we 
will seek to represent the z-transform in terms of 
entities for which we know the z-transform pair. For 
example, we know that the z-transform of anu[n] is 
11—az-—1, |z|>a. So if we can represent some 
signal's z-transform in terms that look like 11— az 


— 1, we will be on the right track to finding the 
inverse. 


Key to this process is the fact that if a z-transform is 
rational (and it will be, for virtually all of our cases 
of interest), then its numerator and denominator can 
be factored according to its poles and zeros. This 
means that an X(z) of the form: 


X(z)=b0+b1z-—-1+b2z-—-2+::+bM 
z-Mlictalz-lt+a2z-—-2+°+aNz— 
N 


can equivalently be expressed as: 


X(z) =zN — M(1 — €1z—1)(1 — G2z—1)::-(1 —CMz-— 1) 
(1 —plz—1)(1— p2z—1)--(1—pNz—1). 


Note how this expression is starting to look similar 
to 11 -—az—1. We can do a little more work on X(z) 
to help with things even more. A rational function 
can be expressed in terms of products of 11—az—1 
terms, but even better, it can also be represented as 
a sum of such terms. Therefore, any X(z) of the form 


X(z)=zN-—-M(1-¢€C1z-1)(1-€2z- 
1) 22 eM 1 Cle pie — 1) po 
ZS A L-S pNoe =k. ) 


can be re-arranged to be expressed as 


X(z)=C0+C11—p1z—1+C21—p2z—1+-- 


+CN1—pNz—1. 


Now, this is something we can work with. Using the 
linearity property of the z-transform, we can take 
the inverse transform of that whole expression to 
yield (assuming the ROC of X(z) extends outwards 
from the outermost pole) 

x[n] = C08[n] + Clp1nu[n] + C2p2nu[n] + --: 
+CNpNnul[n]. 

The plot of x [n] , the inverse z-transform of X ( z 
J=5+2z2-1-2z2-21-162z-1-162z 
= 2. |e |e 2 


Partial Fractions 


Of course, we need to somehow find those C 
constants to be able to know the inverse transform. 
The process of finding X(z)=C0+ C11 —plz 
~1+C21—p2z—1+--+CN1—pNz-1 

from X(z)=b0+b1z—1+b2z—2+::-+bMz 
—M1+alz—1+a2z—2+:::-+aNz—N is known as a 
partial fraction decomposition. 


The first step of the partial fraction decomposition is 
to make sure the order of the denominator 
polynomial is greater than that of the numerator. If 
not, then we first need to carry out a polynomial 
long division and then continue with the remainder. 
For example, suppose 


X(z)=5+2z-—-1-—-z-—-21-16z-1-1 


6Z—-2 , |z|>12 
The long division is a single step: 


61-16z-1- 
16z-2 )5+2z-1-z-2 6-12 
SS 2 Nh ae Sd 


So, 


X(z) =5+2z-—1-—z-—21-—16z-—1-—16z-—2=6+ 
—1+3z-—11-—16z-—1-—16z-2. 


We can now work on that second term, as the 
denominator's order is greater than the numerator. 
The next step, then, will be to factor the 
denominator: 


=l-32=11—102—=1—162—=2=>=1--32 
21001221 )(bei32—1), 


From here, the next step is to express that in terms 
of a sum: 


—1+3z-10—-12z-1)0 +13z-1)=Cl11—-12z 
—1+€C21+13z-—1.The above form assumes that 
each of the poles is unique; if a pole is repeated, e.g. 
1(1 —az—1)2, then that single term in the product 
would correspond to two in the sum, C11 —az 
—1+C1’(1—az—1)2. At this point, we simply solve 
for the C coefficients: 


=) o7—10) = 2 — 1) 132—1) Cl =i 2z 
=F OO 1372— la Po — 1S] C+ 132 
=a COC 127= 1) 137 — l= (Cl C2) +z 
=1(13G1—1T2C2)Cl+C2= 

—1, 18C1 —12C2=3—-C1=3 , C2= —4. 


Given these coefficients, we now have that 
X(z)=6+ 31 —12z-—1-—41+13z-1 , |z| > 12. 


Having split the transform into a sum through 
partial fractions, the inverse transform is now 
straightforward because (due to the linearity of the 
z-transform) we can take the inverse of each piece 
of the sum. So for: 


X(z)=6+ 31 —12z-—1-—41+13z-1 , |z| >12, the 
inverse transform of 1 is 5[n], the inverse of 
11—12z-1 is (12)nu[n], and the inverse of 
11+13z-1 is (—13)nu[n]. For the latter two terms, 
the ROC of the transform is must-know information; 
if the ROC extended inwards instead of outwards, 
the inverses would have been left-sided instead of 
right-sided. Putting all of this together through 
linearity, we have that the inverse transform is: 


x[n] = 66[n] + 3(12)nu[n] —4(—13)nu[n]. 


x{n| 


Discrete-Time Filtering and Filter Design 


Having considered discrete-time signal types, 
systems and their properties, and various signal 
transformations that can assist with analysis, we are 
ready to put this knowledge to use. After all, signals 
and systems are not just interesting mathematical 
entities--they are actually useful! Signals represent 
real-world information, and systems can perform a 
variety of important modifications on them. 

The filtering operation of discrete-time LTI systems 
is similar to the frequency tuning of musical signals 
with a graphic equalizer. [https://www.flickr.com/ 
photos/12577725@N05/6103827986] Yamaha 
EQ-500 by Touho_T, used under [http:// 
creativecommons.org/licenses/by/2.0/] CC BY / 
cropped from original. 


Filtering: LTI Systems in the Frequency 
Domain 


Before we can design a filter to meet some kind of 
desired specification, we first need to remind 
ourselves of what LTI systems can do. We will 
explain filtering from the standpoint of infinite- 
length signals, but the principles apply similarly to 
finite-length signals as well. Recall these two 
important facts: first, that discrete-time signals can 
be represented as an integral of complex sinusoids, 
and second, that complex sinusoids are eigenvectors 


of LTI systems. Putting those facts together, we see 
that an LTI system can be understood as taking in 
discrete-time signals as inputs, and then modifying 
them by scaling their frequency components. The 
amount the system scales each frequency 
component w is the system's frequency response 
H(@). The frequency response is itself simply the 
DTFT of the system's impulse response h[n]. 


So, LTI systems work by taking in input signals and 
scaling each input frequency according to the 
system's frequency response in order to produce the 
output. In this way, they can be thought of as being 
(infinitesimally) finely scaled versions of a graphic 
equalizer, which modifies input music signals 
according to frequency bands. 


Since there is an infinite number of discrete-time 
frequencies @ (as it is a real number between — st 
and st), and each one can be scaled in an infinite 
number of ways (since the scaling factor is any 
complex number), there are obviously an infinite 
number of potential filters/systems. That said, filters 
are usually categorized as being one of four 


different classes. 

The frequency response magnitude and impulse 
response of low-pass filter. Note how the lowest 
frequencies (those near to wo =0) are scaled by a 
(relatively) constant value, while the highest 
frequencies (those near w= +z) are attenuated. A 
high-pass filter does the opposite of a low-pass 
filter; it attenuates lower frequencies and allows 
higher frequencies to pass. Plotted above are the 
frequency response magnitude and impulse response 
of a band-pass filter. These filters allow a narrow 
(contiguous) band of frequencies to pass (here, those 
near w= +72) while attenuating the rest. A band- 
stop filter, like the band-pass filter, also focuses on 
a narrow band of frequencies. But in this case, it 
attenuates that band, while allowing other 
frequencies to pass. 


Filter Types 


The most common filter types are named according 

to the signal frequency components they either stop 

(scale by zero) or allow to pass. For the simplest 

manifestations of these signal classes, the 

magnitudes of their frequency response and their 

corresponding impulse responses are plotted below. 
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As can be seen in the frequency responses above, 
the names of the typical filter classes give a 
straightforward description of the filters' function. 
Note, though, that there is some level of ambiguity 
among the filter names. For example, a low-pass 
filter is also a type of band-pass filter (with low 
frequencies being the "band" that passes), and for 
that matter is also a band-stop filter (the high 
frequencies being the "band" that is stopped). We 
could say the same thing about high-pass filters, as 
well. In practice, then, band-pass and band-stop 
filters typically are understood to pass or stop bands 
that are not about 0 or x, but some other frequency. 
Additionally, the pass- and stop-bands of these two 
filters are typically much narrower than those for 
high-pass or low-pass filters. 

The frequency response magnitude and impulse 
response of an ideal low-pass filter. Note how the 
lower frequencies (those whose magnitude are less 
than the "cutoff frequency" of wc) are scaled by a 
constant value of 1, while the highest frequencies 
(those with magnitudes above the cutoff frequency) 
are completely attenuated. The performance 
constraints for filters are typically expressed as 
allowable deviations from |H(o)|=1 in the filter's 


pass-band (which is indicated by wp) and from | 
H(@)|=0 in the stop-band (the "cutoff" of which is 
indicated by we). 


Filter Design 


Note the simplicity of the four filters plotted above; 
their impulse responses have only 2 or 4 non-zero 
values, and even those are either 1 or -1. In practice, 
discrete-time filters can have impulse response 
lengths that are much longer and much more 
complex. The advantages of such filters is that their 
frequency pass- and stop-points are much sharper, 
with more level pass-bands. We can see this trade- 
off when comparing the frequency and impulse 
response of the low-pass filter above with the 
following "ideal" low-pass filter: 
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We can see why this filter is called "ideal." For all of 
the frequencies below a certain point we, the filter 
scales them by 1, meaning it does not change their 


magnitude at all; then, for all frequencies higher 
than wc, it perfectly zeroes them out. Of course, this 
comes at a cost: the impulse response of this filter 
(h[n] = 2csin(wen)ocn) is infinitely long and is not 
causal (meaning the system needs to know future 
values of the input to calculate the current value of 
the output). Furthermore, because its transfer 
function is not a rational function of finite order, it 
has infinite complexity, and thus cannot be perfectly 
implemented in practice. 


The task of filter design is to work with these two 
constraints, having a desirable frequency response 
without having an overly complex impulse response 
(which affects the complexity and cost of the filter's 
real-world implementation). The "desirability" of a 
filter's frequency response is typically expressed by 
giving certain constraints on its performance: 


|H(w)| 
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Within given filter design parameters, the simplest 
possible implementation is sought. There are two 
main approaches towards implementation, IIR and 
FIR; we will subsequently consider the 


characteristics, advantages, and disadvantages of 
each. 


IIR Filter Design 
Implementation of a discrete-time LTI system. 


Review of IIR Systems 


As discrete-time LTI filters are simply discrete-time 
LTI systems, they can be broadly classified--like 
systems--as being either FIR or IIR. 
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The figure above shows the implementation of any 
discrete-time LTI system; the way all such systems 
vary are simply in their order (the size of N and M) 
and in their coefficients. Systems are IIR if they 
have any non-zero a coefficients. If we see how such 
systems are represented in z-domain, 


H(z)=Y(z)X(z)=b0+b1z2-—-1+bD2z 
—-2+°:-bMz—-Mlt+alz—1+a2z-—-2+ 
“+aNz—-N=2z2N-—-M(z-C1)(z-€2) 
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then we can also say that IIR systems are those 
which have (non-zero and non-infinite) poles. 


IIR Filters 


In discrete-time, the design and implementation of 
IIR filters is largely a matter of appropriating tried 
and tested design elements from continuous-time 
filters. Signal processing obviously predates the 
advent of discrete-time signals and systems, and 
before this time efforts were made to improve the 
design of analog filters. Theory eventually 
approached a point where high-performing filters 
could be created through a straightforward 
mathematical formulation. For discrete-time 
purposes, these filters can be simply translated from 
the Laplace (continuous-time) to the z domain 
through the use of some kind of transformation 
mapping (such as the bilinear transform, where 
s=cz—1z+1). Among the various analog filter 
types that were invented and then eventually 
employed for discrete-time use, we will consider 
low-pass filter implementations of three of the most 
significant types: Butterworth, Chebyshev, and 
Elliptical. 

The impulse response, frequency response 
magnitude, and pole-zero plot of an order N=6 low- 
pass Butterworth filter. 


Butterworth Filter 


The Butterworth filter's defining characteristic is its 


flatness in both its pass- and stop-bands. 
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As with the other IIR filter implementations that 
will be considered here, the filter above has an 
order of just N=6: it takes merely 6 poles and 6 
zeros to create this particular filter, which is 10-100 
times smaller an order than that of similar FIR 
filters. More than other IIR filters, the Butterworth 
also achieves a near-uniform flatness in its pass- and 
stop-bands. However, this does come at a cost; 
compared to other IIR filters of the same order, the 
Butterworth's transition region (centered about w= 
+312 above) is more gradual, meaning there is a 
larger range of frequencies that are not either 
completely passed or attenuated out. 

The impulse response, frequency response 
magnitude, and pole-zero plot of an order N=6 
Elliptic filter. 


Elliptic Filter 


The Elliptic filter has characteristics that are 
essentially the opposite of that of the Butterworth 


filter. 
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Unlike the Butterworth filter, the magnitude of the 
Elliptic filter's frequency response has considerable 
variation about its targets of |H(@)|=1 in the pass- 
band and |H(w)|=0 in the stop-band. However, in 
return for having these ripples across the response, 
note how it makes up for the Butterworth's main 
disadvantage. In contrast to the Butterworth filter's 
somewhat lazy slide from the pass-through region to 
the attenuation region, the Elliptic filter's frequency 


response has a very steep transition at its cutoff 
point; virtually all frequencies reside in either the 
filter's pass- or stop-band. 

The impulse response, frequency response 
magnitude, and pole-zero plot of an order N=6 low- 
pass Chebyshev filter. The impulse response, 
frequency response magnitude, and pole-zero plot of 
an order N=6 Chebyshev Type-1 filter. 


Chebyshev Filters 


The Butterworth filter has a smooth frequency 
response magnitude, but a wide transition region 
between its pass-band and stop-band. The Elliptic 
filter has a very narrow transition region, but must 
trade-off for this with ripples across its frequency 
response. We might wonder if there is a filter that 
could be a kind of compromise between the two, not 
too ripply, and with a decently steep transition at 


the cutoff point. It just so happens that there is. 
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The Chebyshev Type-1 filter has some ripple, but 
note that is only in the pass-band, whereas it has a 
very smooth response at about |H()|=0 in the 
stop-band. Additionally, it has a fairly steep 
transition between the pass- and stop-bands. There 
is a second type of Chebyshev filter that has the 
same steep transition region, while swapping which 


region has ripples and which one is flat: 
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A superposition of frequency response magnitude of 
the four IIR filters considered: Butterworth (black), 
Elliptic (green), Chebyshev Type-1 (blue), and 
Chebyshev Type-2 (red). The phase of the four 
different IIR filter frequency responses: Butterworth 
(black), Elliptic (green), Chebyshev Type-1 (blue), 
and Chebyshev Type-2 (red). The non-linear phase 
response of Elliptic filters is evident when 
comparing the above input and output. 


IIR Filters: An Evaluation 


Among IIR filters, we have seen that they all have 
their own advantages and disadvantages. For a 
given filter-length, the different implementations 
trade-off on having either a smooth frequency 
response magnitude or sharp transition region 
(compromise between the two). 
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As mentioned above, all of these filters also have an 
advantage when compared to FIR filters, for they 
achieve desired frequency response characteristics 
while being relatively small in length (thus less 
hardware/processing is required for their 
implementation). But as you might expect, they do 
have a couple downsides. First: unlike FIR filters, 
IIR filters have poles. The sensitivity of such poles is 
what allows the filters to have such strong effects 
with small filter orders, but that also means that 
deviations in pole locations (which might happen in 
a real-world setting) introduce more distortion than 
a similar deviation of a zero. And second, IIR filters 
have non-linear phase responses. 


Consider especially the pass-band of the filters 
above. While the Butterworth and Chebyshev 
Type-1 filters come closer than the others to being 
linear, all are definitely "curvy" (i.e., nonlinear) 
across the pass-band. A linear phase response is 
desirable because it means that all input frequencies 
will have the same time delay through the filter; a 
non-linear response means that some frequencies 
will be output before others. 


Suppose we have a signal composed of two 
sinusoids, 

x[n] = sin(27501000n) + .25sin(2711501000n), and 
we run this signal through the low-pass Elliptic filter 
shown above (in green). Both of these sinusoids are 
within the pass-band of the filter (as their 
frequencies are +110x and +310zt), so we would 
expect the output of the filter, y[n], to be essentially 
unchanged, perhaps only a bit attenuated because of 
the ripple in the frequency response. However, this 
expectation does not take into account the phase 
response. Because of the phase non-linearity, the 
two sinusoids are shifted different amounts of time, 
resulting in them being out of sync when compared 
to the original plot: 
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FIR Filter Design 
Implementation of a discrete-time LTI system. 
Implementation of an FIR discrete-time LTI system. 


FIR Systems 


Recall the system representation of discrete-time LTI 
systems: 
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FIR systems are those that do not have poles/ 
feedback, so the representation can be simplified: 


From the straightforward time domain expression of 
this kind of system (y[n] = bOx[n] + b1x[n 


—1]+b2x[n-2]+-:-bMx[n—M]) we have the 
following transfer function: 


H(z) = Y(z)X(z) = b0 + b1z—1+b2z—2+:--bMz 
—M=z-—M(z-(C1)(z— (2)::-(z-— CM). 


In the z-domain, such systems have only zeros (not 
counting poles at z=0, or in the case of acausal 
systems, poles at z=). 

FIR filters are designed to have certain tolerances in 
the pass-band and stop-band. 


FIR Filters 


As FIR filters are FIR discrete-time LTI systems, their 
design is simply a matter of adjusting the b 
coefficients of the time-domain representation (or 
equivalently, the location of the zeros in the z- 
domain). As with IIR systems, these systems are 
designed to perform to particular specifications in 
the frequency domain: 


|H(w)| 


FIR filters must be of much higher order (i.e., have 
more non-zero coefficients in their time-domain 
representations) than IIR filters to achieve similar 
performance, but in exchange for this cost have 
several benefits over IIR filters. First, an FIR filter is 
guaranteed to be BIBO stable. Recalling that an LTI 
system is BIBO stable if and only if its impulse 
response is absolutely summable, the finite length of 
FIR impulse responses guarantees this condition. 
And second, unlike the non-linearity with IIR filters, 
FIR filters can be easily designed to have 

a generalized linear phase response (a phase 
response that is linear, with one caveat; see below). 


In order to have a generalized linear phase response, 
FIR filters of length N must be designed so that their 
impulse responses h[n] are either even or odd 
symmetric about the center time value(s) ((N—1)/2 
for odd n, n=N/2 and n=N/2-1 for even N). This 
design can be attained for N that are either even or 
odd. We will consider when N is odd, i.e. when 
N=M-+1 for some even number M. For such a filter 
length, evenness implies that h[n] =h[M-—n]. 
Suppose the impulse response is designed so that 
such is the case. Let's take a look at the system's 
frequency response, using that even symmetry to 
simplify things: 


H(@) = Xn = OMh[nJe —jon=xXn=0M/2-1h[n]e 
—jon+h[M/2]e—joM/2+ Xn=M/2+1Mh[nje 
—jon=Xn=0M/2-—1h[nle—jon+h[M/2]e 


—jaoM/2 + %n=M/2+ 1Mh[M-—nl]e 
—jon=Yn=0M/2-—1h[nJe—jon+h[M/2]e 
—jaM/2 + Xr=0M/2-— 1h[r]e —jo(M—r)=h[M/2]e 
—jaM/2 + %Xn=0M/2—1h[n](e —jon+ ejao(n 
—M))=(h[M/2] + Xn = 0M/2 —-12h[n]cos(w(n 
—M/2)))e—j@M/2=A(o)e—jaM/2. 


Take a close look at the A(@) term. It is called the 
filter's amplitude, because its amplitude is the same 
as the amplitude of the frequency response: |A(o)| 
=|H()|. It is also purely real-valued. Thus the 
frequency response of the filter is a real number 
multiplied by an exponential function e —jawM/2 
that has linear phase. This means that the frequency 
response has a perfectly linear phase, except when 
the sign of A(w) changes. However, since A(@) will 
hover around the value of 1 in the filter's pass-band, 
the filter will have linear phase in the pass-band. 
FIR filter design amounts to creating an A(@) to fit 
particular specifications. 


Creating FIR Filters 
For an FIR filter of odd-length (N=M-+1, where M 
is even) and even symmetry about n= M2 


(h[n] =h[M—n]) we have that the filter's frequency 
response is: 


H(o) = A(w)e—jaM/2, 


where 


A(@)=h[M/2] + *n=0M/2-—12h[n]cos(w(n 
—M/2)). 


Because (generalized) linear phase is desirable, we 
can focus on meeting filter specifications in terms of 
modifying A(o): 


The goal will be to meet the requirements with as 
small of a filter size as possible. To that end, there 
are two pieces of information that will be very 
useful for us. Both involve the "ripples" of A(w) 
centered around A(@)=1 in the pass-band and 
A(@)=0 in the stop-band. It turns out that among 
all filters that meet some defined specifications on 
A(q@) (the error bounds in the pass- and stop-bands 
and the pass and stop frequencies), the filter of 
shortest length M+1 will: 


* be an equiripple filter, meaning that all of the 
ripple oscillations will be of the same deviation 
about A(w)=1 in the pass-band and A(w)=0 in 
the stop-band (in this sense, they "spread" the 
error out equally over the entire frequency 


range) 
* will have A(w) touch the error bounding box 
M2 + 2 times 


Putting those properties together, we can see that 
there are two ways of going about FIR equiripple 
filter design. One way might be to specify the 
desired filter length N=M-+1 (or, 

equivalently, filter order M), and then use the 
properties to create an equiripple filter based upon 
pass-band and stop-band frequencies; because of the 
first property, this filter will minimize the maximum 
deviation from 1 in the pass-band and from 0 in the 
stop-band. To design such a filter, one could use the 
command firpm in MATLAB. Another approach 
could be to specify these desired deviations, and 
then go about finding the filter of minimum length 
that will satisfy these deviations. To do that, one 
could use the firpmord command in MATLAB. 


Equiripple Filter of Length N= 21 

Suppose we wanted to create a low-pass filter of 
length N=21 (order M= 20), with a pass-band 
from @=0 to w =.3x and a stop-band from 

@ = .35n to ®=2. Among all possible FIR filters of 
length N=21, an FIR equiripple filter will have the 
smallest maximum deviation from the optimal 
frequency response magnitudes (of |H()|=1 in 
the pass-band and |H(@)|=0 in the stop-band). 


ccording to the equiripple filter properties, its 
frequency response will touch these error bounding 
boxes M2 + 2=12 times. 

n optimal M= 20 order filter meeting the given 
pass-band and stop-band cutoffs. Note that it 
touches its error bounding box (the red dots on the 


plot) M2+2=12 times. 
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The equiripple filter is also designed, as described 
above, to have an impulse response that has even 
symmetry about its n=10 center time value. As a 
result, it will have a generalized linear phase: 
Equiripple filters can be designed to have impulse 
responses with even (shown here) or odd 


symmetry. The benefit of such design is that the 
corresponding frequency response has generalized 


linear phase. 
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Finally, as this filter is an FIR filter, it will not have 
any poles (except for those at z=0 and, if it were 
acausal, at z= ~). Since its order is M= 20, it does 
have 20 zeros: 

FIR filters only have zeros, and thus are guaranteed 
to be BIBO stable. 


Im(z) 


quiripple Filter of Length N=101 

he filter in the first example, while optimal for 
its length, still has considerable deviation from 
the ideal frequency response of 1 in the pass-band 
and O in the stop-band. Having a more ideal 

esponse is straightforward, simply increase the 
length of the filter; as described above, the 
MATLAB command firpmord could be used to tell 

s the minimum length required to achieve a 
certain margin of tolerance from the ideal 

esponse. Suppose we had some requirement in 

ind, and that command told us a filter length of 
101 was required. The resulting equiripple filter 
of that length is below. 

n FIR filter of length N=101 N=101. Longer 


filters allow for smaller ripples and a sharper 


transition from pass-band to stop-band. 
|H(w)| 


Inverse Filters and Deconvolution 


We have already seen how discrete-time LTI systems 
can be represented mathematically in both the time 
domain (through the impulse response h[n]) and 
frequency domain (through either the frequency 
response H(o)) or the transfer function H(z)). But 
systems can be thought of not only as entities that 
we create to modify signals, but as real-world 
environments which might perform unwanted (or 
only temporarily wanted) changes to the signal. For 
example, a television or cell phone signal traveling 
in the air through a downtown area might reach its 
destination with added echoes. Or a camera lens 
might introduce unintentional blur on a photograph. 
Or a signal might be modified to help with its 
transmission. In each of these cases, we can 
represent the modifications as a system that takes an 
input and produces an output, and we ultimately 
would like to undo the changes the system made; 
we would like to recover the input, given the 
output. A system that seeks to recover this original 
input is called an inverse filter, and the process of 
inverse filtering is also known as deconvolution. 


Inverse Filtering in the z-Domain 


Suppose there is an LTI system G that takes an input 
signal x[n] and produces output y[n]. What we 


would like to do is create another system H that 
takes this y[n] and then, ideally, produce x[n] as its 
output; we will call the actual output x*[n]. 


It might not be clear at first how to create a system 
that recovers x[n] from y[n], especially if we 
consider the problem in the time domain: if 

y[n] =x[n]*g[n], then what must h[n] be so that 
h[n]*y[n] =h[n]*(x[n]*g[n])=x[n]? By the 
associative and commutative properties of 
convolution, you could say we would want 
h[n]*g[n] =5[n] because 6[n]*x[n] =x[n], but even 
then, how do you create such an impulse response? 


The problem is more straightforward when we 
consider it from the z-domain: 
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So if we can create H(z) such that it equals 1G(z), 
then we can recover the original input x[n] by 
filtering the output y[n] through a system whose 
transfer function is H(z). If the original filter G(z) is 
a rational function with certain poles and zeros, 
then the inverse filter will also be a rational 
function; and, since H(z) =1G(z), it follows naturally 
that H(z) will have poles at the locations of G(z)'s 
zeros and zeros at the location of G(z)'s poles: 
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Approximate Inverse Filters 


As remarkable as this simple inverse filtering 
solution is, there is a huge proviso to it: the H(z) 
must be a stable system. If it is also to be causal 
(which is of practical importance), that means that 
all of its poles must lie within the unit circle. As its 
pole locations are simply the zero locations of G(z), 
this means that G(z) is invertible if all of its zeros 
are within the unit circle. If, on the other hand, G(z) 
has zeros on or outside the unit circle, then the 
inverse filter H(z) = 1G(z) cannot be counted on to 
recover the original input. 


There is an intuitive reason why a G(z) with a pole 
on or outside the unit circle could not be inverted. 
Suppose the system G has this simple input/output 
relationship: y[n] =x[n] —x[n—1]. Such a system is 
a primitive high-pass filter, with a single zero on the 
unit circle at z=1. What this zero does is completely 
eliminate the w=0 (i.e., DC) frequency component 
of the input signal. If all we have, then, is y[n], we 
have absolutely no way of knowing what that value 
was, therefore there is no H(z) that could possibly 
recover an arbitrary x[n] from y[n] (of course, if it 
so happened that the original signal x[n] did not 
have any wo =0 component, then it could be 


recovered). Looking at it from a z-domain point of 
view, there would have to be a pole for H(z) at 
z=1, which is BIBO unstable (any input with a non- 
zero @ =0 component will "blow up"). 


Even if G(z) is not invertible, we can still try to 
recover x[n] as best as we can. As we do so, we will 
usually want to make sure our inverse filter H(z) is 
BIBO stable. If it happens that 1G(z) has poles on or 
outside the unit circle, we will simply move them 
inside the unit circle. A straightforward way of 
doing this is by regularization. We choose the 
smallest constant value of r that will make the 
following Ha(z) BIBO stable: Ha(z) =1G(z) +r. 


Matched Filters 


Most discrete-time filters aim to modify an input 
signal according to desired performance in the 
frequency domain; for example, a low-pass filter 
attenuates the high frequency components of an 
input signal. However, there may be instances in 
which the goals of the filter are best understood in 
the time domain. One example of such a filter is the 
matched filter. 


The objective of a matched filter is very 
straightforward: find the location where a particular 
signal occurs within another (larger) signal. While 
this is a task often employed in children's games (I 
Spy) or puzzles (Where's Waldo), it of course has 
more "grown up" uses, as well: for example, sonar 
and radar work by transmitting a signal and then 
listening (searching) for an echoed version of it. 


We have already considered a signal processing 
technique that works very well at finding one signal 
in another. Recall what the inner product operation 
does. The inner product of two signals x and y, (x,y) 
= Ynx[n]y*[n], measures the their correspondence: 
the more alike the two signals are, the greater the 
inner product will be. The more unlike they are, the 
smaller it will be, with a minimum absolute value of 
zero (which means the signals are orthogonal). 


So it makes sense, then, that the process of matched 


filtering will use inner products. But it will take 
more than a single inner product, for we are looking 
for where a small signal occurs in a larger one. What 
we would like to do is perform an inner product 
between the small signal at every possible time 
location in the larger signal, and then examine 
where in the large signal this inner product is at its 
greatest. If we would like to find where the signal 
x[n] occurs in signal y[n], then we would simply 
calculate (y[n],x[n—m]) at every time instance m 
and note where the maximum value occurs. 
Suppose we would like to know where the signal 
x[n] occurs within y[n]. Visually, we can see the 
location is n=15, but to solve the problem 
algorithmically, we would want to take the inner 
product of the small signal against y[n] at every 
time point. 


x{n] 
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n 
y(n] = a[n — 15] + e[n] 


If the process of taking an inner product at every 
time location in a signal sounds familiar, it is 
because that is essentially what discrete-time 
convolution is! Noting that the usual definition 
switches the m and n variables, recall the 
convolution sum: 


y*x = Xny[n]x[m—n]. 


At every time location in the signal, the convolution 
sum is performing an inner-product like operation 
between y[n] and x[—n]. The only thing that is 
missing for it to be an inner product is complex 
conjugation, so technically the convolution sum is 
performing an inner product between y[n] and x*[m 
—nJ; it is determining how similar y[n] and x*[—n] 
are at time m. So if we would like to find where 
x[n] occurs in y[n], then we would simply need to 
carry out the convolution y[n]*x*[—n] and find its 
maximum value. 

The peak value of d[n] =y[n]*x*[—n] indicates 
where the signal x[n] occurs within y[n]. 
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y(n] = a[n ia 15] + e[n] 
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Since we are looking for the maximum value, or (if 
we believe the signal may occur more than once in 
the larger signal) for some threshold to be met, 
matched filtering can be used in the presence of 
noise. In the example above, x[n] appears in a noisy 
signal y[n], but the convolution sum does not 
mind--it simply finds for us the location in y[n] with 
the strongest correlation to x[n]. Because matched 
filtering works in this way, it is also sometimes 
referred to as cross correlation (especially in 
statistical contexts). 


